It should be noted that the Monte-Carlo 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 MC-Potfit 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 Monte-Carlo 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 Monte-Carlo points is not increased appropriately.
A great advantage of using Monte-Carlo 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.
MC-Potfit 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, one mode will be artificially contracted for compatibility
with MCTDH without loss of precision.
a help text is printed:
mcpotfit86 -h
Purpose: creates a natural potential fit.
Usage: mcpotfit<vers><d> [-h|-?] [-ver -rd -w -D name] inpf
vers: program version (e.g. 86 for version 8.6)
d: indicates debug version
inpf: input file (with or without extension ".inp")
optn: -h : print this help-text.
-? : print this help-text.
-ver : Version information.
-rd : Reading of the dvr-file is enforced.
-w : overwrite enabled.
-D name: 'name' denotes the directory where files are
written to, (name in ~.inp file ignored).
-up iup: The integer parameter 'iup' is passed to userpot where
it may be used to select between PES and DMS components.
The order of the options does not matter, but "inpf"
must be the last argument.
xyz-section and ends with
end-xyz-section. At present, the MC-Potfit
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 Monte-Carlo method, No. of points, etc. |
| OPERATOR |
|
Which surface to be used, energy cut-offs. |
| PRIMITIVE-BASIS |
|
Definition of primitive basis. Optional if an external DVR file is specified in the RUN-SECTION. |
| NATPOT-BASIS |
|
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 name-string 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 primitive-basis-section 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, name-directory/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). |
| no-timing |
|
The timing file is not opened. |
| no-OMPpotential |
|
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. |
| no-omega |
|
Do not allocate the SPP-configuration-sampling 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 Ω: (Nsample-coeff × Nconfig) where Nconfig is the product of the number of SPP given in the NATPOT-BASIS-SECTION. Note: Ω is distributed among MPI tasks in case of distributed memory parallelization. Default: not set. | no-omega-t |
|
Relevant only when invert-method = conjgrad is used: then do not allocate the transpose of the SPP-configuration-sampling 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 invert-method = conjgrad is set to reduce cache misses. Also, ΩT is distributed among MPI tasks in case of distributed memory parallelization. Default: not set. |
| no-omega-spp |
|
Do not allocate the grid-sampling matrix ΩSPP.
Per default the grid-sampling matrix ΩSPP
is allocated and the potential is first sampled on all
primitive grid points of a mode for all sampling points
before the integration over the sampling points is
performed. If no-omega-spp is set, the
ΩSPP-matrix is not allocated but 1-D
potential cuts for each sampling point are
calculated and added to the reduced density sequentially.
This saves considerable amounts of memory
but it may also take somewhat longer.
Size of ΩSPP: Ngrid × Nsample, where Ngrid is the number of grid points of a (combined) mode and Nsample is the number of sampling points used to calculate the SPP. Note: ΩSPP is distributed among MPI tasks in case of distributed memory parallelization. Default: no-omega-spp is not set. |
| density = S | If the reduced densities of the modes are calculated (and not read-in, see below) the densities are stored in a file such that they can be re-used 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 read-in. In this case the format is determined from the files directly. | |
| save-evecs | 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 re-used in case one needs to re-calculate the coefficients or add more SPP. | |
| invert-method = 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 which uses a Cholesky decomposition of ΩT×Ω. 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 overlap matrix ΩT×Ω. 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, where R times the trace of ΩT×Ω is used asregularization. The eigenvalues and their regularization are stored for later reference in a separate file. Default is direct. |
|
| cg-tolerance = R | When invert-method = conjgrad is used, cg-tolerance is the length of the residual vector(s) or "error term" divided by the length of B (see invert-method above). Conjugate gradients is considered converged when the length of the residual divided by |B| is less than cg-tolerance. 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 cg-tolerance. Default 1.E-8. | |
| cg-maxiter = I | When invert-method = conjgrad is used, then at most I iterations are used to achieve the requests error-estimate. Default 1000. | cg-precon = S [,I| S] | When invert-method = 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 sub-diagonals used in the band matrix. cg-precon = band,1
will use a tri-diagonal matrix, while cg-precon = 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 comma-separated 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 around 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. |
sampling-only |
|
Do not create a natpot file but only calculate the sampling points and exit. | sampling-energies |
|
Write out the energies that correspond to the sampling points. Default: not set. |
| 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. |
| subtract-pes = S |
|
Subtract a potential stored in a PES file (created with the -pes option of MCTDH)
from the potential that is defined in the OPERATOR-SECTION.
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 cluster-expansion) and only wants to produce a natpot of the difference to the exact potential. If subtract-pes is given, then the OPERATOR-SECTION itself may not use the "pes-file" option but only "usersurf" or built-in surfaces from the MCTDH operator library. subtract-pes cannot be used with extended grids. |
| auto-extgrid |
|
If symmetry operations are present (see SYMMETRY-SECTION below) that do not stay on the DVR grid but lead to grid points outside the primitive grid defined in the PRIMITIVE-BASIS-SECTION, then if auto-extgrid is set, extended mode-grids will be automatically generated such that the symmetry operations stay on an extended primitive grid and a usual grid-based 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. |
| density=S |
|
Save the density in format S, where S either ASCII or binary. Default: ASCII |
| binary-samplings |
|
Save the DVR index points and energies in binary form. Default: ASCII |
| ascii-samplings |
|
Save the DVR index points and energies in ASCII form. This is default. |
Options in the Sampling-Section 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 . |
Monte-Carlo samplings are used for tree different purposes or tasks in MC-Potfit:
SAMPLING-SECTION
sampling-spp = <method 1>
sampling-spp = <method 2>
...
sampling-coeff = <method n>
sampling-coeff = <method n+1>
...
sampling-test = <method N>
end-sampling-section
This will create several sub-sets of sampling points (see examples below) and subsequently
combine them into one large set to perform the calculations. This can for instance be used
to add small fractions of uniform sampling to a metropolis sampling or to combine sets of
metropolis samplings with different temperatures.
In addition, if there are two or more test-sets, detailed statistics
will be written in a separate 'natpot-statistics' file for each subset in the
test-trajectory.If the calculation is a new run, either sampling or sampling-spp and sampling-coeff must be given. The test sampling is optional and only used to estimate error measures of the potfit. A test is in any case performed with the sampling used for calculating the coefficients. A test with an independent set of sampling points is, however, highly recommended. Comparison of a test-set with the coefficient-set is a good measure for the convergence in terms of numbers of sampling points (if the sampling methods are the same).
If the run is a re-run of a previous calculation one can skip the SPP sampling. In this case the SPP are read-in from evec-files if they exist, or are calculated from previously stored density matrices. If sampling-spp is anyhow specified, the density matrices are re-calculated. 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 NATPOT-BASIS-SECTION 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:
SAMPLING-SECTION
sampling-coeff = <method 1>
sampling-coeff = <method 2>
...
sampling-spp = sampling-coeff
sampling-test = <method N>
end-sampling-section
The following examples demonstrate the use of the
sampling[-S1] keyword.
Example 1
SAMPLING-SECTION
sampling-spp = uniform, 10000
sampling-coeff = metropolis, 10000, 1000, 10, 400.0, cm-1
sampling-test = readidx{/path/to/indexfile}, 15000
end-sampling-section
Example 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.
SAMPLING-SECTION
sampling = uniform, 10000
end-sampling-section
Example 3
Metropolis sampling for the coefficients plus an additional small fraction of uniform sampling. For the SPP the same sampling points as for the coefficients are used. The test is done with an independent metropolis sampling.
SAMPLING-SECTION
sampling-coeff = metropolis, 100000, 1000, 10, 400.0, cm-1
sampling-coeff = uniform, 1000
sampling-spp = sampling-coeff
sampling-test = metropolis, 100000, 1000, 10, 400.0, cm-1
end-sampling-section
| 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 = pes-file{path/to/pes} in which case an existing operator file is re-fitted. 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=rd-tfac*rv
r2=rv
3D:
r1=sqrt(rd**2+(tfac*rv)**2
-2d0*tfac*rd*rv*ctheta)
r2=rv
r3=sqrt(rd**2+((1d0-tfac)*rv)**2
+2d0*(1d0-tfac)*rd*rv*cos(theta))
For a homo-nuclear 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 cut-off for exact potential energy surface. All potential energy values greater than R are set to R. |
| vcut > R |
|
Energy cut-off 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 Primitive-Basis-Section.
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 cannot be chosen freely. It must always contain the E operation, the identity. The first column E serves as a reference for all further symmetries.
The columns of the symmetry table under the other symmetry labels implement the symmetry operations. In the column under the label E one must arrange all DOF labels from the Primitive-Basis-Section. The order can be different from the PBASIS-SECTION. If the order of DOFs in a second column is different from the order in the E-column, this is interpreted as an exchange of coordinates within this symmetry operation.
If symmetries are specified one should carefully check the sym.log file. It contains detailed information on the internal representation and mappings of the symmetry-operations. Also a number of tests are performed and logged there to detect if the symmetries are implemented correctly and if the final natpot file is symmetrized correctly.
SYMMETRY-SECTION E SXY1 -------------- x y y x end-symmetry-sectionThe symmetry operation 'SXY' interchanges the grid points of DOF 'x' and 'y'. This exchange does not adhere to any coordinate ranges. If the coordinate 'x' is has grid points in the interval, say, [-1,1], and 'y' in the interval [3,5], it will interchange the first grid point (x1 = -1) of 'x' with the first grid point (y1 = 3) of 'y' (and vice versa), the second grid point of 'x' with the second grid point of 'y' and so forth.
Optionally one may invert the grids. Using the example above with 'x' being defined in [-1,1], and 'y' in [3,5], The following symmetry operation 'SXY2' will map the first grid point (y1 = 3) of 'y' to the first grid point (x1 = -1) of 'x' but the last grid point (xn = +1) of 'x' to the first grid point (y1 = 3) of 'y' and so forth. That is, the 'x' grid is reversed before mapping to 'y'. This again does not adhere to any coordinate ranges.
SYMMETRY-SECTION E SXY2 -------------- x y y -x end-symmetry-section
{±<DOF-Label>±R}
The opening bracket must be followed by a coordinate label which may optionally bear a sign
(scaling of the coordinate is not supported), followed by a constant shift which may be an
arithmetic expression of real numbers.
Besides real numbers also the symbols pi and 2pi can be used.Other than grid-based operations coordinate-based operations do adhere to coordinate values: The grid points of the respective DVR are transformed to coordinate values, altered through the expression in the curly brackets and transformed back to grid points which are ultimately used to perform the symmetrization. Unless one uses extended grids (see below) one therefore needs to make sure that the expression in curly brackets applied to a primitive grid point again leads to a primitive grid point.
The following examples implement the same symmetry operations as SXY1 and SXY2 from the examples above
in coordinate-based form.
SYMMETRY-SECTION
E SXY3 # same as SXY1
--------------
x {y-4}
y {x+4}
end-symmetry-section
SYMMETRY-SECTION
E SXY4 # same as SXY2
--------------
x {y-4}
y {-x+4}
end-symmetry-section
One can define one ore more DOF as periodic with the "periodic" keyword:
SYMMETRY-SECTION
<SYMMETRY-TABLE GOES HERE>
periodic = DOF1 [,DOF2 [,DOF3,...]]
end-symmetry-section
The example below shows an example of the use of the "periodic" keyword:
SYMMETRY-SECTION
--------------------------------------------
E C2 C4 C4^3
--------------------------------------------
x -x y -y
y -y -x x
phi {phi+pi} {phi-pi/2} {phi+pi/2}
--------------------------------------------
periodic = phi
end-symmetry-section
Note: Since periodic DOF typically use exponential or FFT DVRs, mcpotfit will issue
a warning if a DOF with another DVR type is marked as periodic, but will continue anyways and
apply periodic boundary conditions.
The following rules apply for user-implemented symmetries:
SYMMETRY-SECTION
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+2pi} {-la+2pi} {la} {-lb+pi} {lb+pi} {lb+pi} {-lb+pi}
lb {-lb} {lb} {-lb} {la-pi} {-la+pi} {la-pi} {-la+pi}
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
end-symmetry-section
To be able to still perform the symmetry operations one can use the keyword auto-extgrid in the RUN-SECTION. If auto-extgrid is set and it is detected that any symmetry operarion leads to points outside the original primitive grid, these extra points will be automatically added to the respective mode-grids.
When extended grids are used, it is important that all symmetry operations of a point group are implemented, such that that symmetry operations are closed. Otherwise the forward-backward mapping of extended grid points would be incomplete and symmetry operations applied twice would again lead to points outside the grid.
Extended grids can be used for all modes.
compile [options] mcpotfit
compile -P mcpotfit
this will result in the executable mcpotfit86P
compile -m mcpotfit
The resulting executable will be mcpotfit86.mpi-compiler.
MPI support can be combined with OpenMP, of course.
compile -S -m mcpotfit
Note that ScaLAPACK requires MPI. If the -S flag is used, the executable will be
named with the 'S' suffix, for instance
mcpotfit86S.mpi-gfortran
With the default setup 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 older systems, especially if ScaLAPACK is not self-compiled 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.