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 MonteCarlo 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 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 readin and used for generating cpd files.
MCCDP is developed along the lines of MCPotfit 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 nonparallel
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 helptext. ? : print this helptext. c : continue (optimization). mnd : Make name directory. 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). 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 OpenMPparallel 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.
xyzsection
and ends with
endxyzsection
. At present, the MCCPD
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. 
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 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). 
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. 
invertmethod = 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 CPDfit. 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, E_{n} → E_{n} + R tr{A} *exp(E_{n}/R tr{A}) is used as regularization. S=eigen is numerically more costly, but usually (not always!) yields better results. 

loadcuts (= S)  Load the PEScuts that have been created in an earlier run. Note that the cuts must match the samplings defined in the samplingsection. It is recommended to load the samplings from file as well. (See readidx on SamplingSection). The optional S can be a directory in which the program will look for cutfiles. If ommited, S will be the namedirectory.  
tryreusecuts  Try reusing existing PEScuts in namedirectory that have been created in an earlier run. This is useful if a run crashed during cutcalculations and some files have already been created. Note that the cuts must match the samplings defined in the samplingsection. It is recommanded to load the samplings from file as well.  
samplingonly 

Do not create a cpd file but only calculate the sampling points and then exit. 
cutsonly 

Do not create a cpd file but only calculate the PEScuts from a given sampling and then exit. 
samplingenergies 

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. 
autoextgrid 

If symmetry operations are present (see SYMMETRYSECTION below) that do not stay on the DVR grid but 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. 
initguess = S(,S1) 

if S is set to "random", the initial guess basis functions are created from a randomnumber generator, if S = "cuts", the basis functions are created from the 1mode 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 "overrideinitrank" can be used to override the initial rank of "growrank" option. In this case the rank of the cpd file will be the start rank. Default: "random". 
saveinitguess 

The initial guess that is fed into the ALS algorithm is saved on file. Note, that setting saveinitguess 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 growrank (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 growrank is used, 'rank = I' must be set. 
growrank = 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 growrank option. Where: I1: initial rank of the small CPD that is optimized, I2: deltarank, 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 'growrank' is used. Added terms with the 'growrank' option will be initialized with zero coefficients but random basis functions. 
finalizecoeff = 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 finalizecoeff one can
specify if and how the optimal coefficients shall be computed after the
ALS has been completed. If finalizecoeff = no is given, finalization is not performed and the coefficients are used as they come out of the ALS routine. If finalizecoeff = montecarlo is set, the coeffsampling given in the SAMPLINGSECTION is used to recalculate the coefficients allone (this is default). finalizecoeff = cuts is set, the coefficients are recalculated on all cuts that have been loaded or computed. This requires samplingcoeff = samplingspp in the SAMPLINGSECTION.> 
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 . 
MonteCarlo samplings are used for tree different purposes or tasks in MCCPD:
SAMPLINGSECTION samplingspp = <method 1> samplingspp = <method 2> ... samplingcoeff = <method n> samplingcoeff = <method n+1> ... samplingtest = <method N> endsamplingsectionThis will create several subsets 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 testsets, detailed statistics will be written in a separate 'cpdstatistics' file for each subset in the testtrajectory.
If the calculation is a new run, either sampling or samplingspp and samplingcoeff 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 testset with the generatingset is a good measure for the convergence in terms of numbers of sampling points (if the sampling methods are the same).
Note: if PESCuts from a previous calulation are used via the 'loadcuts' keaword, then the SPPsampling 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:
SAMPLINGSECTION samplingcoeff = <method 1> samplingcoeff = <method 2> ... samplingspp = samplingcoeff samplingtest = <method N> endsamplingsectionThe 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
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.
SAMPLINGSECTION samplingcoeff = metropolis, 100000, 1000, 10, 400.0, cm1 samplingcoeff = uniform, 1000 samplingspp = samplingcoeff samplingtest = metropolis, 100000, 1000, 10, 400.0, cm1 endsamplingsection
Example:
MODESSECTION Q1, Q2 # mode 1 Q3 # mode 2 Q5, Q4 # mode 3 Q6 # mode 4 endmodessection
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 nonswapping mode. Therefore one nonswaping 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 errorextimates
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 MODESSECTION 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 = 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 (,S) 

Energy cutoff 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 cutoff 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 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 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 PrimitiveBasisSection. The order can be different from the PBASISSECTION. If the order of DOFs in a second column is different from the order in the Ecolumn, 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 symmetryoperations. 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.
SYMMETRYSECTION E SXY1  x y y x endsymmetrysectionThe 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 (x_{1} = 1) of 'x' with the first grid point (y_{1} = 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 (y_{1} = 3) of 'y' to the first grid point (x_{1} = 1) of 'x' but the last grid point (x_{n} = +1) of 'x' to the first grid point (y_{1} = 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.
SYMMETRYSECTION E SXY2  x y y x endsymmetrysection
{±<DOFLabel>±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 gridbased operations coordinatebased 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 coordinatebased form.
SYMMETRYSECTION E SXY3 # same as SXY1  x {y4} y {x+4} endsymmetrysection
SYMMETRYSECTION E SXY4 # same as SXY2  x {y4} y {x+4} endsymmetrysection
One can define one ore more DOF as periodic with the "periodic" keyword:
SYMMETRYSECTION <SYMMETRYTABLE GOES HERE> periodic = DOF1 [,DOF2 [,DOF3,...]] endsymmetrysectionThe example below shows an example of the use of the "periodic" keyword:
SYMMETRYSECTION  E C2 C4 C4^3  x x y y y y x x phi {phi+pi} {phipi/2} {phi+pi/2}  periodic = phi endsymmetrysectionNote: 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 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+2pi} {la+2pi} {la} {lb+pi} {lb+pi} {lb+pi} {lb+pi} lb {lb} {lb} {lb} {lapi} {la+pi} {lapi} {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 endsymmetrysection
To be able to still perform the symmetry operations one can use the keyword autoextgrid in the RUNSECTION. If autoextgrid 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 modegrids.
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 forwardbackward 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 mccpdthis will result in the executable mccpd86P
compile m mccpdThe resulting executable will be mcpotfit86.mpicompiler. MPI support can be combined with OpenMP, of course.
compile S m mccpdNote that ScaLAPACK requires MPI. If the S flag is used, the executable will be named with the 'S' suffix, for instance
mccpd86S.mpigfortranWith 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 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.
We highly recommend linking sophisticated Lapack packages when building mccpd