Projection Documentation

General Remarks

projection is used to project out some degrees of freedom of the PES (by multiplying the PES with a projector function and integrating over the DOF) to yield a PES of lower dimensionality. This "projected potential" must then be cast into natural potential form in a subsequent potfit run to be usable in mctdh. This feature can be used e.g. to create a multipole expansion of the angular part of a PES, or to create an effective potential by averaging over some DOFs.

Generally, you will use the following procedure:

  1. Write an input file for the projection run.
    Run projection.
    This will yield several vpot_LABEL files, containg the projected potentials.
  2. For each LABEL:
    1. Write an input file to cast the corresponding projected potential into a natural potential representation (normal potfit run). Compared to the input file for the projection run, you must omit from the PRIMITIVE-BASIS-SECTION all the DOFs that were projected out. In the OPERATOR-SECTION, use "pes=none".
    2. In the RUN-SECTION, use "readvpot=/path/to/vpot_LABEL" to specify the vpot file which was created in the projection run. (Alternatively, copy the vpot_LABEL file to the name-directory for the new potfit run, and rename it to just vpot.)
    3. Run potfit.
  3. When writing the operator file for mctdh, in the lines where the projected potentials are referenced, use (for each projected DOF) the operator function that corresponds to the DOF's projector function.
    See the example below.

Command Line Options

projection accepts the following options:

    -h   : print this help-text.
    -?   : print this help-text.
    -ver : Version information.
    -rd  : Reading of the dvr-file is enforced.
    -gv  : force generation of full-grid vpot file.
    -w   : overwrite enabled.
    -mnd : create name directory if it does not exist.
    -D   : 'name' denotes the directory where files are written to,
           (name in ~.inp file ignored).

Input File

The input file for projection follows the same division into sections as the potfit input file. There are only slight differences.

Sections

The various sections are listed below. Those with a STATUS of C are compulsory, O marks optional sections.

XXX Status Description
RUN
C
What files are to be opened, etc.
PROJECTION
C
Definition of projectors for projection runs.
OPERATOR
C
Which surface to be used, coordinate systems, energy cut-offs, etc.
PRIMITIVE-BASIS
C
Definition of primitive basis.

Below, tables of keywords are given for some input section.

A first table describes 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 = S,R' indicates that the keyword takes two arguments, the first is a string, the second a real number. The status column indicates whether the keyword is compulsory, C, or optional, O.

A second table defines the default values for optional keywords and arguments.


Run-Section

Keyword Status Description
name = S
C
The output files will be written to the directory S. Note: If the -D name option is used (see command line options), the name-string given in the input file is ignored.
title = S
O
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)
O
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
O
The DVR file will be generated unconditionally. (This is the default)
deldvr
O
The DVR file will be deleted at the end of the calculation.
readvpot (= S)
O
The vpot file S will be read. S will never be overwritten or deleted irrespective of all other keywords. "readvpot" differs from "readdvr": First, S is a FILE, not a directory. Second, the potential parameters written at the head of file S do NOT overwrite the parameters found in the .inp file. Instead, both sets are compared to each other. If they are not identical, the program stops with an error message. The same happens if the vpot cannot be read. If S is not given, name-directory/vpot is assumed..
genvpot
O
The vpot file will be generated unconditionally. A vpot file specified via the "readvpot" keyword, however, will never be overwritten.
delvpot
O
The vpot file will be deleted at the end of the calculation unless it has not been specified via the "readvpot" keyword.
overwrite
O
Any files already in name directory may be overwritten. (It is safer not to use overwrite but the option -w).
output = S
O
The output will be written to the file name/output. If this keyword is omitted, output will be to the screen. The string S may take the values short or long. When long is specified, additional information might be printed. short is default, i.e. output is equivalent to output=short .
Note: output is default.
screen
O
The output will be directed to the screen, no output-file is opened. Alternatively to screen one may give the keywords no-output or nooutput.
timing
O
Program timing information will be written to the file name/timing   (default).
no-timing
O
The timing file is not opened.

For more details on the format of the various data files see the Potfit Output Documentation.

A projection run allocates a single double precision vector v of full grid size which holds the original potential on the original grid, and which is written to a file "vpot". Additionally, a double precision vector vproj is allocated for the projected potentials. Its size is that of the largest remaining grid (i.e. the grid where all the projected DOFs are omitted). If the projection option error is specified, as many vproj vectors need to be stored as projectors were defined. Note that projection runs cannot make use of "single-precision".


Projection-Section

For the PROJECTION-SECTION, these are the keywords which control the calculation and output of additional files besides the vpot_LABEL files:

Keyword Status Description
error [ = S ]
O
An error measure will be calculated and written to the file S (the default is "projerr.vpot"). This measure is calculated as follows:
Let p be the index of the projector, x be the remaining DOFs, y be the projected DOFs, fp be the projector function, and gp be the complementary projector function (i.e. such that
  ∫dy fp(y) gp(y) = 1 ).
In this notation, the projected potentials are
  Vp(x) = ∫dy V(x,y) fp(y) .
The error measure then is
  ΔV(x) = { ∫dy | V(x,y) − ∑p Vp(x) gp(y) |2 / ∫dy }1/2 .
Of course, this only works if all projectors project onto the same set of DOFs (which is not necessary otherwise). ΔV(x) will be stored as a vpot file and can then be plotted with the showpot utility.
ascii
O
The projected potentials will also be written to ASCII files, named "srf_LABEL.asc". The files contain one energy per line, simply looping over the remaining grid (first remaining DOF runs fastest). This format is usable as input to MCTDH's readsrf statement.

Input syntax

The definition of the projectors themselves is done with the following syntax:

  PROJECTOR LABEL
    MODELABEL  function  parameters
    ...
    option = value
    ...
  end-projector

where:

You can project onto several degrees of freedom at once, by giving multiple "MODELABEL function parameters" lines in the PROJECTOR statement.

List of possible options:

Option Description
scale = R multiplies the resulting projected potential with the overall scale factor R (a real number).

Available functions

In the following table, x denotes the coordinate of the projected DOF. The integration is carried out by a simple quadrature, using the grid points that are specified by the DOF definition in the PRIMITIVE-BASIS-SECTION. This method will yield best results if you use the "recommended DVR" which is also specified in the table. The normalization factors are chosen such that you don't normally need any further factors in the operator file.

Function Parameters Description recomm. DVR Remarks
leg l  m Plm(cos x)
Normalization factor: (2l+1)/2 * (l-m)!/(l+m)!
associated Legendre polynomial of cos(x)		 Leg in .op file, use aslegth:l_m
cos m cos(mx)
Normalization factor: 1/2π
 exp in .op file, use cos[m,0]
gauss a  d exp(−((x−a)/2d)2)
Normalization factor: 1/(sqrt(2π)d)
Gaussian with center a and width d.		   in .op file, use gauss[1/4d2,a]
read filename Read function from the ASCII file filename.
The file must contain one function value per line and have exactly as many lines as the number of grid points for this DOF. It is assumed that the function itself is normalized.		    
How to add new projector functions

After the necessary function parameters, it is possible to specify additional options which modify the projection behaviour:

Option Description
weights_included (Can be abbreviated to wi.)
During integration, the "raw" function values are normally multiplied by the DVR weights. This is simply to take into account the differential volume element. However, if for some reason you know that the raw function values already include the DVR weight, you must specify this option to prevent the normal inclusion of the weights during integration.
Example: the veigen keyword in MCTDH can be used to write eigenfunctions of 1D operators to a file. These eigenfunctions already include the square root (!) of the DVR weights. So if you used the square of such an eigenfunction as input, you would need wi.
root_of_weights_included (Can be abbreviated to rwi.)
Like "weights_included", this tells the integration routine that the function values already include the square root of the DVR weights.
no_weights_included (Can be abbreviated to nwi.)
This option is just available for completeness and tells the integration routine that the read function values are "raw", i.e. they do not include the DVR weights. This is also the default.


Operator-Section

Keyword Status Description
pes = S {pesopts}
C
S is the name of a potential energy surface which must have been encoded in the mctdh package. 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/Liouvillian Documentation -- Available Surfaces. If you want to fit a vpot file that was created by a different program (e.g. adproj, or projection itself) use the "pes = none" option. This means that potfit uses only the vpot file and does not look for a potential energy surface.

Example:
For the lsth surface use 'pes = lsth {tfac = 0.5d0}'. tfac is a mass-dependent parameter which is needed to transform between binding and Jacobian coordinates. Let r1,r2, and r3 denote the binding and rd,rv, and theta the Jacobian coordinates. Then the transformation rules are as follows:

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 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 centre 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.
order = S/I
O
order defines the ordering of the variables passed to the potential routine. If S=dof, then the ordering is as in the PRIMITIVE-BASIS-SECTION. Note that this is the default, order=dof need not to be given. If S=mode (or S=particle), then the ordering is as in the NATPOT-BASIS-SECTION. Finally, one may give a blank separated list of integers, which define the ordering with respect to the dof-ordering.
E.g.:     order = 3 1 2
passes the third dof as the first argument to the potential routine. The first dof becomes the second and the second dof the third argument. The arguments of order must be a permutation of the numbers 1...ndof. See the log-file to inspect if the order is correct.
vcut < R
O
Energy cut-off for exact potential energy surface. All potential energy values greater than R are set to R.
vcut > R
O
Energy cut-off for exact potential energy surface. All potential energy values less than R are set to R.
cspot = R1,R2,R3
O/C
Whether an additional centrifugal potential of the following form is added to the potential energy surface (so far only installed for the 3D LSTH-surface in Jacobian coordinates):
cspot(rd)=[R1(R1+1)-2*R2^2]/[2*R3*rd^2]

where R1 denotes jtot, R2 denotes K and R3 denotes mass_rd. 'cspot' is compulsory if option '-jtot' has been specified on the command line.


Primitive-Basis-Section. Ordering of arguments

The primitive basis is defined in the same way as it is done in the input file of an mctdh run (see the MCTDH Primitive-Basis-Section ). There is, however, one important difference. The order of the degrees of freedom defines the ordering of the arguments passed to the potential routine (unless you use the order keyword). To see, how a particular pes is implemented see the file funcsrf.F (just type   mcb funcsrf.F ).  Inspect also the log-file log, which is created when running potfit. There one finds a message like: 

BKMP2 PES : 3D model in Jacobian Coordinates
Dissociative Coordinate : rd
Vibrational Coordinate  : rv
Angular Coordinate      : theta
A line like Dissociative Coordinate : theta would hint to a wrong ordering of the degrees of freedom in the Primitive-Basis-Section.
The arguments passed to the potential routine may be reorderd with the aid of the order keyword of the Operator-Section.
NB: In MCTDH the ordering of the degrees of freedom in an HAMILTONIAN-SECTION (i.e. the ordering of the mode-line) must be consistent with the potential definition. There the ordering within the Primitive-Basis-Section is arbitrary.


Example Input Files

This example shows how to produce the "BMKPE" expansion of the BMKP surface for H4, as suggested by Pogrebnya & Clary [Chem. Phys. Lett. 363 (2002) 523, Section 2.2].