It should be noted that the MonteCarlo integration error manifests mainly in an integral over the residue potential (the part not represented by an optimal natpot once the SPP are determined). In Particular this means: if the potential can be completely expressed in the SPP basis obtained by MCPotfit in the first step, then also the obtained natpot file will be numerically exact. If this is not the case, that is, if the SPP cannot represent the exact potential, one needs to keep in mind that the role of the SPP is actually to interpolate the potential between the MonteCarlo points which were used to calculate the coefficients. This can lead to a paradox situation that the total error of the fit can increase if more SPP are used and the number of MonteCarlo points is not increased appropriately.
A great advantage of using MonteCarlo techniques is that correlated weights are implemented straight forwardly by means of the sampling density, i.e., the distribution of the sampling points. At present, generating a uniform distribution and a Boltzmann distribution using a Metropolis algorithm are implemented. Furthermore, external samplings can be read in and used for generating natpot files.
MCPotfit does not yet support all features of the original Potfit.
Implemented features are outlined below. Another difference is that the
internal structure of the resulting natpot file is somewhat different from the
original Potfit. A contracted mode is not mandatory. When creating the natpot
file, however, there will be an retrospectively contracted mode for compatibility
reasons.
a help text is printed:
mcpotfit84 h
Purpose: creates a natural potential fit. Usage: mcpotfit<vers><d> [h?] [ver rd w D name] inpf vers: program version (e.g. 84 for version 8.4) d: indicates debug version inpf: input file (with or without extension ".inp") optn: h : print this helptext. ? : print this helptext. ver : Version information. rd : Reading of the dvrfile is enforced. w : overwrite enabled. D name: 'name' denotes the directory where files are written to, (name in ~.inp file ignored). The order of the options does not matter, but "inpf" must be the last argument.
xyzsection
and ends with
endxyzsection
. At present, the MCPotfit
input may contain four to six
sections, as listed below. Those with a STATUS of
C are compulsory, O marks optional sections.
XYZ  Status  Description 

RUN 

Defines general parameters. 
SAMPLING 

Defines the MonteCarlo method, No. of points, etc. 
OPERATOR 

Which surface to be used, energy cutoffs. 
PRIMITIVEBASIS 

Definition of primitive basis. Optional if an external DVR file is specified in the RUNSECTION. 
NATPOTBASIS 

Size of natural potential basis, optional contracted mode, etc. 
SYMMETRY 

Define symmetries of the PES fit 
Below, tables of keywords are given for the input sections.
The following tables describe the keywords. The number and type of
arguments is specified. The type is S for a character string, R
for a real number, and I for an integer. For instance,
'keyword = I,R[,S]'
indicates that the keyword
takes two or optionally
three arguments where the optional argument is indicated
by square brackets. Here, the first argument is an integer,
the second a real number, and the optional third one
a character string. The
status column indicates whether the keyword is compulsory, C, or
optional, O.
Keyword  Status  Description 

name = S 

The output files will be written to the directory S.
Note: If the D name option is used (see
To Run the
Program), the namestring given in the input file is
ignored. 
title = S 

The string S is taken as the title of the run and printed to the log and output files. Note: everything that follows the equal sign '=' till the end of the line is taken as title! 
readdvr (= S) 

The DVR information will be taken from the DVR file in directory S. Note that the primitivebasissection of the input file is completely ignored if this keyword is given. A DVR file residing in directory S will NEVER be overwritten or deleted, even if the keywords deldvr and/or gendvr have been specified. If S is not given, namedirectory/dvr is assumed. 
gendvr 

The DVR file will be generated unconditionally. (This is the default) 
deldvr 

The DVR file will be deleted at the end of the calculation. 
overwrite 

Any files already in name directory may be overwritten. (It is safer not to use overwrite but the option w). 
output = S 

The output will be written to the file
name/output . The string S may take the values
short or long. When long is specified,
additional information on grids and weights are printed.
short is default, i.e. output is equivalent
to output=short .Note: output is default. 
dvronly 

The program generates the dvr file and then stops. 
timing 

Program timing information will be written to the file
name/timing (default). 
notiming 

The timing file is not opened. 
noOMPpotential 

If openMP parallelization is used with this flag no calls to the PES routine are performed in parallel. This is useful if the PES routine is not thread safe. 
noomega 

Do not allocate the SPPconfigurationsampling matrix Ω which is needed to calculate the coefficients, but calculate it on the fly as needed. This saves considerable amounts of memory but also takes considerably longer, especially for large sampling sets. Dimensions of Ω: (N_{samplecoeff} × N_{config}) where N_{config} is the product of the number of SPP given in the NATPOTBASISSECTION. Note: Ω is distributed among MPI tasks in case of distributed memory parallelization. Default: not set. 
noomegat 

Relevant only when invertmethod = conjgrad is used: then do not allocate the transpose of the SPPconfigurationsampling matrix Ω^{T}. If Ω is allocated, use Ω instead, otherwise calculate it on the fly as needed. This saves considerable amounts of memory but also takes considerably longer, especially for large sampling sets. Note: Ω^{T} is only allocated and used when invertmethod = conjgrad is set to reduce cache misses. Also, Ω^{T} is distributed among MPI tasks in case of distributed memory parallelization. Default: not set. 
density = S  If the reduced densities of the modes are calculated (and not readin, see below) the densities are stored in a file such that they can be reused in a later run if for instance more coefficients are needed. The density keyword describes the format of the file. S can be either ASCII or binary, where binary is default. The ASCII format is useful if a visualization of the density, for instance with gnuplot is desired. The density keyword is ignored if the reduced densities are readin. In this case the format is determined from the files directly.  
saveevecs  In Addition to the reduced densities (see density above) save a copy of all single particle potentials as well and in the same format as the densities. The SPP can be reused in case one needs to recalculate the coefficients or add more SPP.  
invertmethod = S (,R)  Method used to solve the system AX=B where A is the overlap matrix
(Ω^{T}×Ω), and B is (Ω^{T} V). X are the coefficients
of the SPP in the potfit. S can be one of direct, conjgrad preudo or eigen. When direct is selected, the overlap matrix (Ω^{T}Ω) es explicitly calculated and inverted using the LAPACK DPOSV (or ScaLAPACK PDPOSV) routine. This can be quite expensive but is numerically exact. If conjgrad is selected, a conjugate gradients method is used that avoids explicit calculation of the overlaps. If eigen is selected, the matrix (Ω^{T}×Ω) is diagonalized and inverted using LAPACK DSYEV (or ScaLAPACK PDSYEV) This allows for regularization, with the parameter R, there R times the trace of Ω^{T}×Ω is used asregularization. Default is direct. 

epsinvert = R  Regularization parameter added to the diagonal elements of the SPP overlap matrix. Can be helpful to speed up conjugate gradients.  
cgtolerance = R  When invertmethod = conjgrad is used, cgtolerance is the length of the residual vector(s) or "error term" divided by the length of B (see invertmethod above). Conjugate gradients is considered converged when the length of the residual divided by B is less than cgtolerance. If a contracted mode is present, B has multiple vectors, hence, conjugate gradients is considered converged when the length of the largest residual divided by B is less than cgtolerance. Default 1.E8.  
cgmaxiter = I  When invertmethod = conjgrad is used, then at most I iterations are used to achieve the requests errorestimate. Default 1000.  
cgprecon = S [,I S]  When invertmethod = conjgrad use the preconditioner S, where
S can be one of "diag" in which case a Jacobi preconditioner is used, or "band"
in which case a banded matrix is used or "sparse", see below. S="band"
requires an integer I
specifying the number of subdiagonals used in the band matrix. cgprecon = band,1
will use a tridiagonal matrix, while cgprecon = band,0 is equivalent to the
Jacobi preconditioner.
Note: The behavior of the band matrix preconditioner is somewhat unpredictable as it is not
guarantied that the condition number of Ω^{T}Ω is reduced. It can also get worse.
The band matrix preconditioner can therefore speed
up the convergence substantially, but also reduce it substantially. For metropolis sampling we
always observed some speed up using the Jacobi preconditioner. For uniform sampling, a
preconditioner is usually not needed as Ω^{T}&Omega is close to
unity and well conditioned. The sparse preconditioner uses a sparse matrix as a preconditionier which is computed from selected matrix elements of the SPP overlap matrix. One specifies the number of matrix elements mode wise as a commaseparated list in the form <mode_1>:<N_1>,<mode_2>:<N_2>, etc., where <mode_i> is the mode number and <N_i> is the number of matrix elements for this mode. Note: there must not be spaces arount the ':'. If all <N_i> are equal to the numbers of SPP in their mode, then the preconditioner matrix is not sparse any more but complete and the preconditioner is exact. 

samplingonly 

Do not create a natpot file but only calculate the sampling points and exit. 
iseed = I 

Calculates the initial random seed from I. If not set, the random seed is calculated from the current system time. In case of MPI, the rank number is added in each task and the system time used is the one from rank 0. 
subtractpes = S 

Subtract a potential stored in a PES file (created with the pes option of MCTDH)
from the potential that is defined in the OPERATORSECTION.
S is the full path to a PES file. This option is useful if one already has an approximate description of the potential (for instance a clusterexpansion) and only wants to produce a natpot of the difference to the exact potential. If subtractpes is given, then the OPERATORSECTION itself may not use the "pesfile" option but only "usersurf" or builtin surfaces from the MCTDH operator library. 
autoextgrid 

If user implemented symmetry operations are present (see SYMMETRYSECTION below) and if these symmetry operations lead to grid points outside the primitive grid defined in the PRIMITIVEBASISSECTION, then if autoextgrid is set, extended modegrids will be automatically generated such that the symmetry operations stay on an extended primitive grid and a usual gridbased symmetrization can be performed. The additional grid points are subsequently removed such that the original grid is restored. If not set, an error message is produced if any symmetry operation leeds to outside the primitive grid. Default: not set. 
Options in the SamplingSection are the following
sampling[S1] = S2, I1 [,I2, I3, R [,S3]] 

The sampling method and the number of sampling points to use. See remark below . 
samesets 

Use the same set of sampling points for all tasks. remark below . 
MonteCarlo samplings are used for tree different purposes in MCPotfit:
If the run is a rerun of a previous calculation one can also skip the SPP sampling. In this case the SPP readin from evecfiles if they exist, or are calculated from previously stored density matrices. If samplingspp is anyhow specified, the density matrices are recalculated. Note, that if the densities are read in, at present only the grid dimensions are checked. It is not checked weather the density represents a valid representation of the given mode. That is, interchanging modes in the NATPOTBASISSECTION will not lead to an error if the combined grids of the modes have the same sizes.
The first argument of the sampling[S1] keyword is the sampling method, followed by a set of parameters. Possible sampling methods  with parameters  are:
The following examples demonstrate the use of the sampling[S1] keyword.
Example 1
SAMPLINGSECTION samplingspp = uniform, 10000 samplingcoeff = metropolis, 10000, 1000, 10, 400.0, cm1 samplingtest = readidx{/path/to/indexfile}, 15000 endsamplingsectionExample 2
Uniform sampling with 10000 points for generating the SPP, the coefficients and also for testing the natpot. The sampling points will be different for the three individual tasks, only the method of their generation is the same.
SAMPLINGSECTION sampling = uniform, 10000 endsamplingsectionExample 3
Uniform sampling with 10000 points for generating the SPP, the coefficients and also for testing the natpot. The sampling points will be the same for all three individual tasks.
SAMPLINGSECTION sampling = uniform, 10000 samesets endsamplingsection
Keyword  Status  Description 

pes = S {pesopts} 

S can be either the keyword "usersurf" in which the program assumes
that the potential energy routine has been implemented in the file
$MCTDH_DIR/source/potfit/userpot.f90
or it can be the name of a potential energy surface which has
been encoded in the MCTDH package. A third option is pes = pesfile{path/to/pes} in which case an existing operator file is refitted. The first option has been provided to allow easy implementation of user provided routines using a FORTRAN 90 interface. Note, that options from the input file are not supported in this case, i.e, there must not be any {pesopts}. In case that a potential energy from the MCTDH library the name is interpreted in a case sensitive manner. If a surface depends on additional parameters, one can change the encoded default values by adding those parameters via the string '{pesopts}'. Parameter names specified within the braces '{}' are processed in a case sensitive manner. The selected surface determines which parameter names are possible. For a description of the available surfaces see Hamiltonian Documentation  Available Surfaces. Note that fitting 'vpot' files is not supported at present. Example: 2D: r1=rdtfac*rv r2=rv 3D: r1=sqrt(rd**2+(tfac*rv)**2 2d0*tfac*rd*rv*ctheta) r2=rv r3=sqrt(rd**2+((1d0tfac)*rv)**2 +2d0*(1d0tfac)*rd*rv*cos(theta)) For a homonuclear diatomic molecule tfac=0.5d0, in general tfac=m1/(m1+m2). rv denotes the distance between the two atoms of the diatom, rd the distance between the third atom and the center of mass of the diatom, and theta the angle between rd and rv. Note: tfac is not used if binding coordinates are selected for a calculation. 
vcut < R 

Energy cutoff for exact potential energy surface. All potential energy values greater than R are set to R. 
vcut > R 

Energy cutoff for exact potential energy surface. All potential energy values less than R are set to R. 
BKMP2 PES : 3D model in Jacobian Coordinates Dissociative Coordinate : rd Vibrational Coordinate : rv Angular Coordinate : thetaA line like Dissociative Coordinate : theta would hint to a wrong ordering of the degrees of freedom in the PrimitiveBasisSection.
The symmetries are given as a table of symmetry operations where each column of the table implements one symmetry operation. The first line of the symmetry table contains the labels by which the symmetry operations are addressed. The label can be chosen freely but must be a single word. The only exception is the first label of the symmetry table which must always contain the E operation, the identity. It serves as a reference for all further symmetries.
The columns of the symmetry table under each symmetry label then implement the symmetry operations. The first symmetry operation is always E, the identity. In the column under the label E one must arrange all DOF labels from the PrimitiveBasisSection. The order can be different from the PBASISSECTION, however, all further symmetry operations are implemented with respect to the column containing the E operation. If the order of DOFs in a second column is different from the order in the E column, this is interpreted as exchanges of coordinates within this symmetry operation. Note, that this is done on a gridlevel and not a coordinate level. The coordinate values are not exchanged, but just the grid points. The symmetry operations different from E may not only interchange the order of the DOF but also reverse the the order of the primitive basis points indicated by a "" sign.
SYMMETRYSECTION E SXY1  x y y x endsymmetrysectionto preserve this symmetry in the fit. Note the order of x and y in the second column with respect to column one. Since this operation exchanged grid points, one must take care that the primitive grits are chosen accordingly.
If the PES has the properties V(x,y) = V(y,x) then one might implement
SYMMETRYSECTION E SXY2  x y y x endsymmetrysectionThis operation interchanges the DOF x and y, and also reverses the primitive basis points of the x coordinate.
The following operations are allowed:
The following constraints to the table apply:
The following roles apply for userimplemented symmetries:
SYMMETRYSECTION E C2 SigmaA SigmaB S4 S4^3 C2A C2B  z z z z z z z z R R R R R R R R x usersym usersym usersym usersym usersym usersym usersym y usersym usersym usersym usersym usersym usersym usersym a usersym usersym usersym usersym usersym usersym usersym la la la la lb lb lb lb lb lb lb lb la la la la ua ua ua ua ub ub ub ub ub ub ub ub ua ua ua ua r1a r1a r1a r1a r1b r1b r1b r1b r2a r2a r2a r2a r2b r2b r2b r2b va va va va vb vb vb vb r1b r1b r1b r1b r1a r1a r1a r1a r2b r2b r2b r2b r2a r2a r2a r2a vb vb vb vb va va va va endsymmetrysection
compile [options] mcpotfit
compile P mcpotfitthis will result in the executable mcpotfit84P
compile m mcpotfitThe resulting executable will be mcpotfit84.mpicompiler. MPI support can be combined with OpenMP, of course.
compile S m mcpotfitNote that ScaLAPACK requires MPI. If the S flag is used, the executable will be named with the 'S' suffix, for instance
mcpotfit84S.mpigfortranAt compile time the linker will look in the standard system search path for the file 'libscalapack.a'. If the library is not in the search path or has another name, one can alter this accordingly in the variable MCTDH_SCALAPACK_LIBS of compiler in use in the compile.cnf file. Note that ScaLAPACK uses BLACS for communication, which is usually shipped with the Scalapack package and compiled into the ScaLapack library. On some systems, especially if ScaLAPACK is not selfcompiled but installed from the systems package manager, BLACS may be installed in a separate library, usually called libblacs.a. If this is the case one also needs to add the BLACS library in the compile.cnf file.
At present, ScaLAPACK is only used for solving for the natpot coefficients.