The program starts with an initial guess of the CPD and improves it iteratively using the alternating least squares (ALS) method. As the traditional ALS requires multible integrals over the complete primitive grid, these integrals are replaced by Monte-Carlo integrals. This on the one hand enables the use of larger potentials, but on the other hand also introduces a statistical error into the resulting potential representation. However, the error can be estimated again by statistical methods.
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 cpd files.
MC-CDP is developed along the lines of MC-Potfit and requires a very similar input file.
Furthermore, the program makes heavy use of BLAS/LAPACK libraries
and therefore profits a lot from linking to an optimized BLAS/LAPACK
library.
The use of BLAS/LAPACK libraries can be reduced by issuing the
-noblas command line flag. This could lead to some
performance improvement when linked to a non-parallel
standard BLAS/LAPACK library while using OpenMP at the same time.
It usually makes no sence to use -noblas in a pure MPI
or a serial run.
a help text is printed:
mccpd86 -h
Purpose: creates a canonical poliadic fit.
Usage: mccpd [-h|-?] [-ver -rd -w -c -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.
-c : continue (optimization).
-mnd : Make name directory.
-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.
-noblas : When compiled with -P option (OpenMP parallelization)
this program relies heavily on being linked to an
optimized and OpenMP-parallel BLAS/LAPACK installation.
Use this flag if you are using OpenMP parallelization
but a serial (!) BLAS/LAPACK library is linked.
(The default BLAS/LAPACK library is serial.)
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-CPD
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. |
| MODES |
|
Definition of mode combinations . |
| 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). |
| 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. |
| invert-method = S (,R) | Method used to solve the system AX=B where A is the overlap matrix
of the CPD terms, and B is the overlap of the potential with a CPD term.
X are the coefficients
of the CPD-fit. S can be one of direct or eigen. When direct is selected, the overlap matrix is inverted using the LAPACK DPOSV (or ScaLAPACK PDPOSV) routine which uses a Cholesky decomposition. This is default. If eigen is selected, the matrix is diagonalized and inverted using LAPACK DSYEV (or ScaLAPACK PDSYEV) Note that regularization (R) is essential when solving the linear system to prevent diverging CPD terms. Per default R is set to be the square root of machine precision. If S=direct is chosen, The diagonals of the matrix A are multiplied with (1+R). If S=eigen is chosen, En → En + R tr{A} *exp(En/R tr{A}) is used as regularization. S=eigen is numerically more costly, but usually (not always!) yields better results. |
|
| load-cuts (= S) | Load the PES-cuts that have been created in an earlier run. Note that the cuts must match the samplings defined in the sampling-section. It is recommended to load the samplings from file as well. (See readidx on Sampling-Section). The optional S can be a directory in which the program will look for cut-files. If ommited, S will be the name-directory. | |
| try-reuse-cuts | Try reusing existing PES-cuts in name-directory that have been created in an earlier run. This is useful if a run crashed during cut-calculations and some files have already been created. Note that the cuts must match the samplings defined in the sampling-section. It is recommanded to load the samplings from file as well. | |
| sampling-only |
|
Do not create a cpd file but only calculate the sampling points and then exit. | cuts-only |
|
Do not create a cpd file but only calculate the PES-cuts from a given sampling and then exit. |
| sampling-energies |
|
Write out the energies that correspond to the sampling points. This default for Metropolis sampling. |
| 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. |
| 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. |
| initguess = S(,S1) |
|
if S is set to "random", the initial guess basis functions are created from a random-number generator, if S = "cuts", the basis functions are created from the 1-mode cuts through the PES. Otherwise S is interpreted as a path to a cpd file that is used as initial guess. If S is intepreted as a path to a cpd file, the keyword "override-initrank" can be used to override the initial rank of "grow-rank" option. In this case the rank of the cpd file will be the start rank. Default: "random". |
| save-initguess |
|
The initial guess that is fed into the ALS algorithm is saved on file. Note, that setting save-initguess will trigger the calculation of the CPD coefficients and hence will be numerically approximately as costly as one ALS optimization for one mode. If not set, the CPD coefficients are initialized with zero. They do not enter the numerics. |
| maxiter = I |
|
I is the maximum number of iterations the ALS will perform. Default: I = 1000. |
| rank = I |
|
I is the maximum rank of the CPD. If grow-rank (see below) is not used, I is also the initial and final rank of the CPD. 'rank = I' can be ommited if the initial guess is read from file or for continuation runs. If grow-rank is used, 'rank = I' must be set. |
| grow-rank = I1, I2, I3 |
|
When optimizing a CPD it proved very benefitial to start optimizing a small CPD tensor first
and sequentally growing it, i.e., to repeatedly add terms after a number of ALS iterations.
This can be done with the grow-rank option. Where: I1: initial rank of the small CPD that is optimized, I2: delta-rank, i.e. number of terms terms added after I3 iterations of the ALS algorithm. Adding additional I2 terms is done every I3 iterations of the ALS algorithm until the total rank is equal to the value given in the 'rank = I' option above (or until the number of ALS iterations exceeds the 'maxiter' option). 'rank = I' is therefore mandatory if 'grow-rank' is used. Added terms with the 'grow-rank' option will be initialized with zero coefficients but random basis functions. |
| finalize-coeff = S |
|
The coefficients that come out of the ALS iterations are
not the optimal ones. They minize the CPD on the SPP sampling-
points where one mode is replaced by unity.
With the keyword finalize-coeff one can
specify if and how the optimal coefficients shall be computed after the
ALS has been completed. If finalize-coeff = no is given, finalization is not performed and the coefficients are used as they come out of the ALS routine. If finalize-coeff = monte-carlo is set, the coeff-sampling given in the SAMPLING-SECTION is used to re-calculate the coefficients allone (this is default). finalize-coeff = cuts is set, the coefficients are re-calculated on all cuts that have been loaded or computed. This requires sampling-coeff = sampling-spp in the SAMPLING-SECTION.--> |
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-CPD:
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 'cpd-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 fit. A test is in any case performed with the sampling used for calculating the coefficients or SPP. A test with an independent set of sampling points is, however, highly recommended. Comparison of a test-set with the generating-set is a good measure for the convergence in terms of numbers of sampling points (if the sampling methods are the same).
Note: if PES-Cuts from a previous calulation are used via the 'load-cuts' keaword, then the SPP-sampling must match the sampling the cuts have been created with.
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
Example:
MODES-SECTION
Q1, Q2 # mode 1
Q3 # mode 2
Q5, Q4 # mode 3
Q6 # mode 4
end-modes-section
The internal numbering of the modes is simply done according to the line in which the
mode is specified and is mostly arbitrary. There is one important exception:
If symmetry operations are used and modes are swapped (for instance mode 1 and 3 could
swap upon some symmetry operation in the example above)
then the last mode must be one that only maps onto itself upon all symmetry
operations. It must not be one that swaps.
The reason is, that during the ALS iterations the coeffictions of the CPD are not
obtained correctly when modes are swapped. This is, however, automatically 'healed' when the
ALS algorithm optimizes the next non-swapping mode. Therefore one non-swaping mode
must always be the last one.
(If there are only swapping modes in the system, then
there is no choice and the optimization will be performed anyways but the error-extimates
in the 'iterations' file will be inaccurate. As the coefficients do not enter the
working equations this only affects the output and not the numerics.
The correct coefficients will - for performance reasons - then only
be calculated once in the very end, after the ALS iterations have been completed.)
The MODES-SECTION is ignored in continuation runs or if the initial guess is loaded from file.
| 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/mcpotfit/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 (,S) |
|
Energy cut-off for exact potential energy surface. All potential energy values greater than R are set to R. S denotes the unit of R, default is au. |
| vcut > R (,S) |
|
Energy cut-off for exact potential energy surface. All potential energy values less than R are set to R. S denotes the unit of R, default is au. |
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, mccpd 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] mccpd
compile -P mccpd
this will result in the executable mccpd86P
compile -m mccpd
The resulting executable will be mcpotfit86.mpi-compiler.
MPI support can be combined with OpenMP, of course.
compile -S -m mccpd
Note that ScaLAPACK requires MPI. If the -S flag is used, the executable will be
named with the 'S' suffix, for instance
mccpd86S.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.
We highly recommend linking sophisticated Lapack packages when building mccpd