OZ is a software for numerically solving the Ornstein-Zernike equation with a closure (MSA, PY, HNC, RY, ZH, and RHNC) for homogeneous systems with particles possessing central forces

Version 4.2.2

Documentation update June 2015

Per Linse

Division of Physical Chemistry

Department of Chemistry

Lund University, P.O. Box 124

S-221 00 Lund, SWEDEN

This document can be downloaded from www.physchem.lu.se/sm

The software OZ may not be redistributed or modified, with the exception of the routines in source file ozuser.F90, without the written consent of the copyright owner. The copyright owner does not take any responsibility for any software in the code or its documentation.

We hope that the software will be useful for you and request a proper reference (OZ, Version 4.2.2, P. Linse, Lund University, Lund, Sweden 2015) when using it.

Preface

This documentation describes the software OZ. After a short introduction, information of how to install the software, to set up an environment for a user, and to run it is provided. Thereafter, an explanation of all input variables divided into subsections is given, and the namelist technique used is explained in an appendix. The remaining chapters deal with user-provided routines and more technical aspects important for a more advanced use of the software.

Only the use of the software is given in the documentation. It is assumed that the user has a general understanding of the relevant theory. Very few references are provided here; however, further references are given in the source code, often at the top of the routines.

The software OZ is written in Fortran 95. It can be run interactively or in a batch environment. The connection between external units and files is made in a command file.

Table of contents

1. Introduction...................................................................................................................................... 5

2. Installation........................................................................................................................................ 6

3. Getting started.................................................................................................................................. 7

4. Input................................................................................................................................................. 8

4.1. Namelists................................................................................................................................... 8

4.2. Input variables........................................................................................................................... 8

4.2.1. nmlSystem.......................................................................................................................... 8

4.2.2. nmlPotential...................................................................................................................... 10

4.2.3. nmlIteration...................................................................................................................... 11

4.2.4. nmlInitGam...................................................................................................................... 12

4.2.5. nmlUseBr......................................................................................................................... 13

4.2.6. nmlScattItens.................................................................................................................... 13

5. User-provided routines................................................................................................................... 14

6. Distributed files.............................................................................................................................. 15

7. External units.................................................................................................................................. 16

Routines............................................................................................................................................. 17

7.1. oz.f90...................................................................................................................................... 17

7.2. ozaux.f90................................................................................................................................. 17

7.3. ozlib.f90.................................................................................................................................. 18

7.4. ozuser.f90................................................................................................................................ 18

8. Appendix A. Namelist.................................................................................................................... 19

9. Appendix B. Suggestions or problems?......................................................................................... 20

1. Introduction

OZ is a software for numerically solving the Ornstein-Zernike equation with a closure for homogeneous systems with particles possessing central forces.

An arbitrary number of components can be selected. All potentials that can be described as a sum of A/r**n terms with different prefactors and exponents operating between pairs of components can be used. In addition, hard-sphere repulsion, screened Coulomb potential, and a potential cutoff can be applied. The closures available are the classical (i) Mean Spherical Approximation, (ii) Percus-Yevic, and (iii) Hypernetted chain, some mixed ones as the (iv) Rogers-Young and (v) Zerah-Hansen closures, and (vi) the related Reference Hypernetted chain approach.

The main results provided by the software are (i) thermodynamic properties as internal energy, pressure, and isothermal compressibility, (ii) potential energy, total correlation function, direct correlation function, radial distribution function, and the bridge functions, (iii) Fourier transformation of many structural functions, (iv) scattering intensity, form factors, and total structure factor for multi-shell models, and (v) zero and second moment conditions.

The numerical solution of the Ornstein-Zernike equation with the appropriate closure is obtained iteratively improving the nodal function until self-consistency within a given tolerance is achieved. Two principal solution methods are available; (i) the original Picard method including an extension by Ng and (ii) the method by Labik et al., involving an expansion of the nodal function in a basis set and combined Newton-Raphson and Picard iterations. The long-range Coulomb tails are analytically subtracted before and added after the numerical Fourier transformations. The initial nodal functions can be selected (i) to be zero, (ii) from analytic solutions of the PY approximation of a hard-sphere system, the MSA of a restricted primitive model, or the MSA of a one-component model with a screened Coulomb potential, or (iii) from a previous calculation.

It is easy to add user-provided routines to the software to enhance the capability of it. These routines may, e.g., perform additional analyses or facilitate sequential calculations. The user-provided routines are collected in the file ozuser.F90 and do not belong to the supported software. An existing input variable controls the calls of these routines.

The specification of input variables characterizing the model and controlling the calculation is given in an input file. That makes it easy to make sequential calculations with only minor differences in the specification. Both names and values of the variables are given in input file using the Namelist protocol, making it easy to read and prepare input files. Blocks of input data are collected into separate groups.

Whenever, say two components, species, or states A and B appear ordered as (A, B), the order of the pairs of components, species, or states is (AA, AB, BB), etc. This ordering is important for input variables.

2. Installation

This chapter describes how to install the software OZ on a Unix(Linux)-based workstation after obtaining the installation file oz.xxx.tar.gz, where xxx should be replaced by three digits (version number).

In the following, your top directory is referred to as $(HOME). The installation involves the following steps:

– Create the directory OZ in the directory $(HOME)

– Copy the file oz.xxx.tar.gz to the directory OZ

– Go to file directory OZ

– Type 'tar -zxf oz.xxx.tar.gz'

The tar command (i) creates the subdirectory Src, which in the following is referred to as the source directory, and (ii) places the source files into the source directory. The name and content of the distributed files are given in chapter 6.

– Go to directory Src

– Open the makefile and if necessary edit it for appropriate compiler

– Type 'make'

The make command starts compilation and linking. After completion, (i) additional files have been created in the source directory, (ii) a subdirectory Bin has been created and now contains the executable oz.exe, and (iii) the file oz has been copied from the source directory to $(HOME)/bin.

Note, the environmental variable PATH has to contain $(HOME)/bin for the software to function.

– The installation file oz.xxx.tar.gz in directory OZ can now be removed

The installation of the software is now completed.

3. Getting started

This chapter describes how to execute the software.

– Create a work subdirectory, in which corresponding input and result files will reside

This directory will in the following be referred to as the work directory. The location of the work directory is arbitrary.

– Copy an input file, here referred to as job.in, to the work directory

An input file could, e.g., be retrieved from www.physchem.lu.se/sm. At this moment we will not consider the content and the meaning of the variables in the input file. The input variables are the subject of the next chapter.

– Type 'oz job'

After execution the output file job.out have been created in the work directory.

The content of the output file should be essentially the same as that obtained by execution trough www.physchem.lu.se/sm. If not, probably some error has occurred in the transfer of the code or during the installation.

4. Input

All input data read from the file FIN are done with the namelist procedure. Each namelist contains related input variables. Further information on the namelist procedure is given in Appendix A.

4.1. Namelists

The program uses the following namelists in the order given below.

nmlSystem Main variables

nmlPotential Variables describing the potentials

nmlIteration Variables describing the solution of the equations

nmlInitGam Variables for generation of initial nodal function

nmlUseBr Variables for RHNC calculation

nmlScattItens Variables for calculation of form factor and scattering intensity

4.2. Input variables

This subsection contains all input variables that are used for reading data from file FIN. The variables are grouped together and listed below their namelist name. The first line contains the name, the type, the dimension, and, if any, the default value of the variable. The following lines briefly explain the use of the variable. If the variable only can attain a limited number of values, these are listed. If description about writing data is given, the output is made on file FOUT (if nothing else is stated).

A practical point. Often there is a need of reading a variable conditionally, i.e., if some condition is true the variable is assigned a value and used later. If the condition is false, the presence of the variable in the input file does no harm; thus it has not to be deleted. In the latter case the input file contains redundant data.

4.2.1. nmlSystem

The namelist nmlSystem contains variables that describe the main features of the calculation as the method and closure as well as flags that control output.

txTitle character(80)

User-provided title.

txUser character(10)

Character string for selecting appropriate user-provided routine called from OZUSER. See driver routine OZ from where routine OZUSER is called.

txMethod character(8)

'picard' Iterative solution according to method of Picard

'picardng' Iterative solution according to the method of Picard with Ng's extension.

'labik' Iterative solution according to the method of Labik, Malijevsky, and Vonka.

TxOZEqu character(8) /'analytic'/

'analytic' Analytic solution of the OZ-equation (only for NPT≤2).

'numeric' Numeric solution of the OZ-equation.

txClosure character(4) (1:nptpt)

Character sting for selecting approximation for each pair of components.

'msa' Mean spherical approximation.

'smsa' Soft-sphere mean spherical approximation.

'py' Percus-Yevick approximation g=[exp(-bu)](1+g).

'hnc' Hypernetted chain approximation g=exp(-bu+g).

'ry' Rogers-Young approximation g= exp(-bu){[exp(fg)-1]/f}.

'zh' Zerah-Hansen approximation g={exp[f(-bu+g)]-1}/f+1.

'ref' Reference hypernetted chain approximation g=exp(-bu+g+B R).

f=[1-exp(-r*FAlpha)]**FPower

FAlpha real (1:nptpt)

Factor describing the RY and ZH approximations.

FPower real (1:nptpt) /nptpt*1/

Factor describing the RY and ZH approximations.

NPT integer

No of components.

TXPT character*5 (1:npt)

Labels of components.

Dens real (1:npt)

Number density of the components.

Temp real

System temperature (Kelvin).

nBin integer

No of bins representing the radial part of the real and reciprocal spaces.

NBF integer (1:nptpt)

No of basis functions for each pair of components (only used if txMethod='labik'). Optimal value is 0.05*nBin to 0.1*nBin.

DR real

Grid length of the real space. The grid length of the reciprocal space is p/DR.

iTest integer /0/

0 Nothing. The normal use.

1 Extensive output for maintaining purposes.

iShow integer /0/

0 Nothing.

>0 Calculated functions are listed at every iShow bin.

iPlot integer /0/

0 Nothing.

>0 Calculated functions are plotted.

iList integer /0/

0 Nothing.

>0 Calculated functions are listed at every iList bin on file FLIST.

lSI LOGICAL /.false./

.true. Calculation of form factors and scattering intensity.

.false. Nothing.

lMomCond LOGICAL /.false./

.true. Calculation of zero and second moments (for Coulomb systems).

.false. Nothing.

lHSBr LOGICAL /.false./

.true. Calculation of bridge functions of hard sphere system according to Malijevksy and Labik.

.false. Nothing.

lNumCheck integer /.false./

.true. Calculating and writing

(i) the accuracy of the solution,

(ii) the division of the nodal function in basis functions,

(iii) the subtraction of the long range tail.

.false. Nothing. The normal use.

4.2.2. nmlPotential

The namelist nmlPotential contains variables that describe the potentials between the components.

RCut real /1.0E30/

Cutoff distance of the potential. The cutoff used is min(RCut,DR*(N-1))

DiamHS real (1:nptpt) /3*0.0/

Hard-sphere contact distance.

TAlpha real (1:nptpt)

Factor used to subtract the long-range tail, see IFTSub.

nPot integer (1:nptpt) /3*0/

nPot(IPTJPT) denotes the no of UCoeff(M)/r**iPot(M) terms of component pair IPTJPT.

iPot integer (1:mnPot, 1:nptpt)

iPot(M,I) denotes the exponent of term no M of component pair I.

UCoeff real (1:mnPot, 1:nptpt)

UCoeff(M,I) denotes the coefficient of term M of component pair I. Energy unit assumed is kJ/mol.

lScrC LOGICAL /.false./

Flag for screened Coulomb potential. The 1/r term is multiplied by exp(-r/ScrLen). iPot=1 is required.

.true. Screened Coulomb potential.

.false. No multiplication of the 1/r term.

ScrLen real /0.0/

Screening length of the screened Coulomb potential.

lErfc LOGICAL /.false./

Flag for complementary error function. The 1/r term is multiplied by erfc(r*Alpha); Alpha=UAlpha/RCut. iPot=1 is required.

.true. The complementary error function is used.

.false. No multiplication of the 1/r term.

lErf LOGICAL /.false./

Flag for error function. The 1/r term is multiplied by erf(r*Alpha); Alpha=UAlpha/RCut. iPot=1 is required.

.true. The error function is used.

.false. No multiplication of the 1/r term.

UAlpha real /0.0/

Damping factor used for the complementary error and error functions.

lNeuBackg LOGICAL /.false./

Flag for neutralizing background.

.true. g(r) is replaced with g(r)-1 when calculating the internal energy and pressure.

.false. No replacement.

4.2.3. nmlIteration

The namelist nmlIteration contains variables that control the solution of the equations.

DampPi real /0.0/

Damping factor for the Picard iteration (0=<DampPi=<1). Only used if txMethod='picard' or 'picardng'.

DampNR real /0.0/

Damping factor for the Newton-Raphson iteration (0=<DampNR=<1). Only used if txMethod='labik'.

DampRC real /0.0/

Damping factor for the refinement-cycle (0=<DampRC=<1). Only used if txMethod='labik'.

TolPi real /1.0E-5/

Tolerance of the Picard iteration. Only used if txMethod='picard' or 'picardng'.

TolNR real /1.0E-6/

Tolerance of the Newton-Raphson iteration. Only used if txMethod='labik'.

TolRC real /1.0E-6/

Tolerance of the refinement-cycle. Only used if txMethod='labik'.

nIterPi integer /50/

Maximum no of Picard iterations. If convergence is not achieved a warning message is issued, and the calculation is continued. Only used if txMethod='picard' or 'picardng'

nIterNR integer /20/

Maximum no of Newton-Raphson iterations. If convergence is not achieved a warning message is issued, and the calculation is continued. Only used if txMethod='labik'.

nIterRC integer /15/

Maximum no of refinement-cycle. If convergence is not achieved a warning message is issued, and the calculation is continued. Only used if txMethod='labik'.

IFTSub integer /0/

Flag for selecting function used to subtract the long range tail of the direct correlation function and the nodal function. Only used if Coulomb potential.

0 erfc(r*TAlpha)/r is used.

1 exp(-r*TAlpha)/r is used.

4.2.4. nmlInitGam

The namelist nmlInitGam contains variables that control the generation of the initial nodal function.

txGamInit character(6)

Text label used for selecting procedure of obtaining the initial nodal function.

'OLD' Taken from file FAUX which contains the nodal function from a previous calculation.

'ZERO' Set to zero.

'PYHS' From the solution of the PY approximation of a hard-sphere system.

'MSARPM' From the solution of the MSA of a restricted primitive model.

'MSASCR' From the solution of the MSA of a one-component model with a screened Coulomb potential.

DiamHSInit real /0.0/

Hard-sphere diameter used when txGamInit='phys', 'msarpm', or 'msascr'

lModify LOGICAL /.false./

.true. Modification of the initial nodal function. The modification is controlled by variables ULimit, RLimit, RSwitch, and GamFac.

.false. No modification.

ULimit real (1:nptpt) /(-1.0E+30)/

The nodal function is modified such that g(r)=1 if the potential < ULimit.

RLimit real (1:nptpt) /1.0E+30/

The nodal function is modified such that g(r)=1 if r > RLimit.

RSwitch real (1:nptpt) /1.0E+30/

The nodal function is modified such that g(r)=g(r)(FAUX) if r > RSwitch and equal to the nodal function (read from file FAUX) times GamFac if r < RSwitch.

GamFac real (1:nptpt) /nptpt *1.0/

Factor used when r < RSwitch.

4.2.5. nmlUseBr

The namelist nmlUseBr contains variables that describe the reading and manipulation of the bridge functions for RHNC calculation. The RHNC calculation is only implemented for txMethod='labik'.

FacBeg real /0/

The fraction of the bridge functions used (read from file FBRAUX) in the first refinement-cycle of the RHNC calculation.

FacInc real /1/

The fractional increment of the bridge functions used (read from file FBRAUX) in subsequent refinement-cycles of the RHNC calculation. The total fraction is at most 1.

iSelect integer (1:nptpt) /3*1.0/

This variable makes it possible to perform RHNC calculation for only some component pair(s).

0 Bridge function is zero (HNC).

1 Bridge function according to FacBeg and FACINC.

4.2.6. nmlScattItens

The namelist nmlScattItens contains variables that control the calculation of the form factors and scattering intensity.

nShell integer (1:npt) /npt*1/

No spherical shells with different contrast.

RShell real (1:npt,1:mnShell) /npt*DiamHS/

Radius of the different shells.

CShell real (1:npt,1:mnShell) /1/

Contrasts of the different shells.

5. User-provided routines

The file ozuser.f90 contains routines that are conditionally called and which might be modified by the user. This allows the software to be extended without having to modify the core of it.

The existing code in the distributed file ozuser.f90 reflects the current use of it. The procedure to insert a new routine is as follows:

– At the appropriate position in driver ozuser, extend the test of txUser and write the appropriate call to the new subroutine

– Add your new subroutine in ozuser.f90

– Compile

6. Distributed files

This chapter contains a complete list of the distributed source files.

oz.f90 Main code of OZ

ozlib.f90 Library routines

ozuser.f90 User-provided routines

oz.ver Description of version changes

oz Command file for execution

makefile Makefile for compilation and linking

archive Command file for filing

7. External units

A file is specified by its name and type, which are separated by a dot. OZ uses file type in a predetermined way to describe the use of the file and leaves file name to describe the project.

File type	Generic name	Content
in	FIN	Input data
out	FOUT	Output data
aux	FAUX	Iteration variables
baux	FBRAUX	Bridge functions
list	FLIST	List data
user	FLUSER	Avail for the user

Flow chart of the use of external files. Bold arrows denote files that are always required or generated, while thin arrows denote files that might be required or generated, depending on input variables.

Routines

This chapter lists the routines in the software. The routines are listed according to the Fortran files they are located in and in the order they appear in the files.

7.1. oz.f90

OZModule	main module of OZ
OZ	main routine of OZ
IOSystem	read and write main input data
IOPotential	perform i/o on variables describing the potential and set up potentials
SMSA	calculate parameters for soft-sphere msa
IOIteration	perform i/o on variables describing the solution of the oz equation
IOGamInit	perform i/o on variables setting initial gam and set initial gam
GamAux	it is assumed that the potential is finite (no overflow check)
GetBr	read the bridge functions and gradually increase br
Picard	iterate according to the Picard method
PicardNg	calculate a new gam according to Ng
Analyse	perform the analyses
ScattIntens	calculate the scattering intensity
MomentCondition	reduced zero and second moments as a function of r of a two-comp system
HSBridgeFunc	hard-sphere bridge functions
CheckSolution	check the accuracy of the solution
TestOutput	write test output

7.2. ozaux.f90

Closure	calculate g, h, c, or br from gam by using a closure
dClosure	calculate dc/dgam by using a closure
OZEqu	calculate cft, hft, or gamft from one of the other using the oz equation
OZEquAnalytic	calculate analytically cft, hft, or gamft from one of the other using the oz eq.
OZEquNumeric	calculate numerically cft, hft, or gamft from one of the other using the oz eq.
dOZEqu	calculate dgamft(is)/dcft(it) by using the oz equation
CalcEnergyPressure	calculate potential energy and pressure
EnergyPressureRdf	calculate from (r,g) lists (dist file)
CalcCompressibility	calculate isothermal compressibility from cft
CalcSFPartial	calculate SFPartial from hft
CalcSFTotal	calculate s(tot) from s(partial)
GROptP	print the first maximum and the following minimum of g
GROpt	calculate the first maximum and the following minimum of one g
GRContact	calculate contact value of g
GRFit	calculate extreme values of g and corresponding r
PYHS	solution of a one-component hard-sphere model using the py approximation
MSARPM	solution of a two-component restricted primitive model using the msa
MSAScr	solution of the primitive model with screened Coulomb potential using msa
MSASub	calculate constants and g1 based on eta, gamma, and ak
MSASubFunc	determine roots
RootOk	check if the root is ok, subroutine to MSASub

7.3. ozlib.f90

FFTC	perform a (inverse) ft of a complex function
FFT3D	calculate the 3d ft (real) of a real function using FFTC
FTSub	add h(r) or hft(k)=ft(h(r)) given below to a real function
Erf1	return the error function, erf(x)
Trap	integration by the trapezoidal rule
PolFit	calculate least square fit to a polynomial
PolVal	calculate the value of the polynomial for a given x-value
MatMulx	multiply two matrices, c = a*b
MatInv	matrix inversion with solution of linear equations
Root2	divide an interval into equal segments & search for zero crossings
Root3	find a root of a function known to be there. Bernt's method
VInter	interpolate in a 1d array
LinInter	linear interpolation between (x0,y0) and (x1,y1)
WriteHead	write a heading
WriteDateTime	write date and time
WriteLog	write a logical variable
WriteInt	write a integer
WriteFlt	write a real variable
WriteVec	write a vector
WriteMat	write a matrix
CpuAdd	add and write total cpu time elapsed since start
CpuTot	write total cpu time elapsed since start
Second	return cpu time used since start (in seconds)
Warn	write a warning message
Stop	write a warning message and stop the process
LowerCase	change a string to lower case
UpperCase	change a string to upper case
Plot2	plot a function

7.4. ozuser.f90

OZUser	driver of user-provided routines
LoopOZ	loop over different complete calculations
CompressibilityVirialR	calculate the virial compressibility
UeffMSA	calculate the effective macroion one-component potential (msa approx.)
UeffHNC	calculate the effective macroion one-component potential (hypernetted approximation)
UseUeffHNC	use the effective one-component potential (hnc approximation)
PYHSHW	hard sphere - hard wall pair correlation function (py approximation)
CalcBrSim	calculate bridge functions from simulated rdfs
GDiff	calculate the rms difference between two rdfs
ListRdf	list rdf in appropriate format

8. Appendix A. Namelist

– A namelist consists of (i) one start statement, (ii) one end statement, and (iii) a list (possibly empty) between the start and end statements.

– The start statement consists of an ampersand (&) followed by the name of the namelist.

– The end statement consists of a slash (/).

– The list is composed of variables and their values separated with either by commas.

– The order of the variables and their assignments is normally irrelevant.

– However, if a variable is assigned values more than once, the last one takes place.

– Arrays are allocated either element by element, consecutively by listing values separated by commas, or by a combination of these to ways.

– A value, say 5, appearing r times may be expressed as r*5.

– Note the order Fortran stores arrays elements: the left most index is running fastest and the right most slowest.

– Comments are allowed and should proceed by an exclamation sign (!).

– Generally, if a program reads several namelists, they have to occur in the inputfile in the same order as they are read from the program. Here, this is relaxed; the order of the namelist is arbitrary.

– A namelist with no list (empty namelist) has to be specified.

– Namelists that are not read by the program do no harm.

Example: Consider the following declarations

character(4) :: title

integer(4) :: m, n

real(8) :: ARR1(1:3), ARR2(1:2,1:2), ARR3(1:2,1:2)

The namelist nmlSystem below illustrates an element-by-element array assignment

&nmlSystem

title = 'head',

n = 10, m = 10,

n = 20,

arr1(1) = 1.5, arr1(2) = 2.5, arr1(3) = 2.5,

arr2(1,1) = 1.0, arr2(2,1) = 2.0, arr2(1,2) = 3.0, arr2(2,2) = 4.0,

arr3(1,1) = 1.0, arr3(2,1) = 2.0, arr3(1,2) = 3.0, arr3(2,2) = 4.0,

which could, for example, be shortened to

&nmlSystem

title = 'head',

n = 20,

m = 10,

arr1 = 1.5, 2*2.5,

arr2 = 1.0, 2.0, 3.0, 4.0,

arr3(1,1) = 1.0, 2.0, arr2(1,2) = 3.0, 4.0,

9. Appendix B. Suggestions or problems?

I would be happy to hear from you if you have any suggestions or problems with the OZ software or with this documentation. If so, send

– Date

– Name

– E-mail address

– OZ version number

– Host computer (name and configuration)

– Suggestion or problem in as much detail as possible

to: Per.Linse@physchem.lu.se