PES generators

SURF,Start1D=label,options [surf]

The SURF program allows for the calculation of the potential energy surface around a reference structure as required for the calculation of anharmonic frequencies (see the VSCF and VCI programs). The reference structure is supposed to be a (local) minimum or a transition state of a double-minimum potential. The potential is represented by energy grid points rather than an analytical representation. Within the SURF program the potential energy surface is expanded in terms of normal coordinates, linear combination of normal coordinates or localized normal coordinates. Consequently, a harmonic frequency calculation needs to be performed first. The potential will then be represented by a multi-mode expansion, i.e. a hierarchical scheme given by $$\label{eq:h} V(q_1,\dots,q_{3N-6}) = \sum_i V_i(q_i) + \sum_{i<j} V_{ij}(q_i,q_j) + \sum_{i<j<k} V_{ijk}(q_i,q_j,q_k) + \dots$$ with $$\begin{aligned} \label{diff1} V_i(q_i) & = & V_i^0(q_i) - V(0) \\ V_{ij}(q_i,q_j) & = & V_{ij}^0(q_i,q_j) - \sum_{r\in\{i,j\}} V_r(q_r) - V(0) \\ V_{ijk}(q_i,q_j,q_k) & = & V_{ijk}^0(q_i,q_j,q_k) - \sum_{\stackrel{\scriptstyle r,s\in\{i,j,k\}}{r>s}} V_{rs}(q_r,q_s) - \sum_{r\in\{i,j,k\}} V_r(q_r) - V(0) \qquad \\ V_{ijkl}(q_i,q_j,q_k,q_l) & = & \dots \label{diff4}\end{aligned}$$ where $q_i$ denotes the coordinates. This expansion needs to be terminated after an $n$-body contribution as controlled by the keyword NDIM. The SURF program is fully parallelized in a sense that the calculation of different grid points is send to different processors (embarassingly parallel MPPX scheme). The START1D keyword is mandatory and defines the label where to jump in the input in order to do an electronic structure calculation which is terminated by the SURF command. This way the quality of the potential energy surface is defined. For example, the input for the calculation of a CCSD surface looks like:

label1
hf
ccsd
surf,start1D=label1

The SURF program is based on an iterative algorithm, i.e. grid points will be added automatically to the grid representation of the potential until a convergence threshold will be met. This guarantees a well-balanced description of the different terms in the expansion of the potential and simultaneously minimizes the number of ab initio calculations for a representation of the potential. For further details see:
G. Rauhut, Efficient Calculation of Potential Energy Surfaces for the Generation of Vibrational Wave Functions, J. Chem. Phys. 121, 9313 (2004).
T. Hrenar, H.-J. Werner, G. Rauhut Accurate Calculation of Anharmonic Vibrational Frequencies of Medium Sized Molecules Using Local Coupled Cluster Methods, J. Chem. Phys. 126, 134108 (2007).

The following options are available:

• BATCH3D=n After calculating a number of grid points within the iterative interpolation scheme the convergence of the individual surfaces will be checked and, if provided by the keyword DUMP, dumped to disk. This leads typically to 3-5 iterations and thus the same number of restart points within the calculation of the 1D, 2D, … surfaces. As the number of 3D and 4D terms can be very large this is not sufficient in these cases. Therefore, the lists of 3D and 4D terms is cut into batches which will be processed subsequently. BATCH3D and BATCH4D control the number of 3D and 4D surfaces within each batch. By default BATCH3D is set to 30 times the number of processors and BATCH4D to 10 times the number of processors. Accordingly the number of restart points is increased. Smaller values for BATCH3D and BATCH4D, e.g. BATCH3D=20, increase the number of restart points on cost of the efficiency of the parallelization. Note, this keyword is only relevant for SURF calculations, but not for XSURF runs.
• DELLOG=n For large molecules or in the case of modelling the 3D and 4D terms, the .log-file may become huge. First of all the .log-file can be directed to scratch within the electronic structure program, i.e. logfile, scratch. The option DELLOG=1 always truncates the .log-file in a way that it contains only the very last energy calculation. Default: DELLOG=0.
• EXT12D=value Outer regions of the potential energy surfaces may be determined by extrapolation rather than interpolation schemes. By default extrapolation is switched off, i.e. Ext12D=1.0 and Ext34D=1.0. However, an extrapolation of 10% in case of the 1D and 2D contributions to the potential (Ext12D=0.9) and of 20% in case of the 3D and 4D terms (Ext34D=0.8) may be useful as it usually stabilizes the fitting procedure.
• FIT1D=n The maximum order of the polynomials used for fitting within the iterative interpolation scheme can be controlled by the keywords FIT1D, FIT2D, FIT3D, FIT4D. The default is given by 8. However in certain cases higher values may be necessary, but require an appropriate number of coarse grid points, which can be controlled by MIN1D etc.
• INFO=n INFO=1 provides a list of the values of all relevant program parameters (options). Default: INFO=0.
• MAX1D=n The maximum number of coarse grid points can be controlled by the keywords MAX1D, MAX2D, MAX3D, MAX4D. These 4 keywords determine the maximum number of ab initio calculations in one dimension for each 1D, 2D, 3D and 4D surface. The defaults are currently MAX1D=24, MAX2D=16, MAX3D=10, MAX4D=8. Presently, values larger than 24 are not supported.
• MIN1D=n The minimum number of coarse grid points can be controlled by the keywords MIN1D, MIN2D, MIN3D, MIN4D. These 4 keywords determine the minimum number of ab initio calculations in one dimension for each 1D, 2D, 3D and 4D surface. The defaults are currently MIN1D=4, MIN2D=4, MIN3D=4, MIN4D=4. Presently, values larger than 24 are not supported.
• MPG=n Symmetry of the normal modes is recognized by the program automatically. Only Abelian point groups can be handled at the moment. Symmetry of the modes will be determined even if the NOSYM keyword is used in the electronic structure calculations. In certain cases numerical noise can be very high and thus prohibits a correct determination of the symmetry labels. Symmetry can be switched off by using MPG=1.
• NDIM=n The keyword NDIM=n terminates the expansion of the PES after the $n$-body term. Currently, at most 4-body terms can be included, but the default is set to 3. Please note, when you use NDIM=4 as a keyword for the SURF program, you need to pass this information to the VSCF and VCI programs also. Otherwise these programs will neglect the 4-body terms.
• NGRID=n Based on a coarse grid of ab initio points a fine grid will be generated from automated interpolation techniques. The keyword NGRID=n determines the number of equidistant grid points in one dimension. NGRID=n has to be an even number. The default is currently set to 16. Note that the number of grid points also controls the extension of the $n$-dimensional potential energy surfaces (see keyword SCALE) and thus influences many internal thresholds which are optimized to the default value of NGRID. The number of grid points also determines the number of basis functions in the grid-based VSCF program. At present the maximum grid size is 36.

• ORIENT Allows to specify a certain orientation of the molecule. With ORIENT=yes (Default) the orientation is choosed automatically according to the asymmetric parameter of the molecule. To choose a certain orientation, ORIENT=XC need to be set. X represents a number from 1 to 3 (in arabic or roman letters), and C need to be set to r or l. For example, ORIENT=IIl orientates the molecule according to the IIl convention. ORIENT=old does not rotate the molecule at all.
• PLOT=n PLOT=n plots all nD surfaces and a corresponding Gnuplot script in a separate subdirectory (plots1) in the home-directory in order to allow for visualization of the computed nD surfaces. E.g. the command “gnuplot plotV1D.gnu” in the plots1 directory produces .eps files for all 1D surfaces. Default: PLOT=0.
• SADDLE=n Standard SURF calculations expect the reference structure to be a (local) minimum on the PES, i.e. SADDLE=0 (default). Alternatively, one may start the PES generation from a transition state, which is recommended for the calculation of double-minimum potentials. This situation is not recognized automatically and thus requires the keyword SADDLE=1. Within XSURF calculations, this keyword needs not to be provided.
• SCALE=value The extension of the potential energy surfaces is determined from Gauss-Hermite quadrature points. Using a fine grid NGRID=16 the surface stretches out to the NGRID/2$^{th}$ Gauss-Hermite point, i.e. 4.69, in each direction (see keyword NGRID). As these values are fairly large within the calculation of fundamental modes, a scaling factor, SCALE=f, has been introduced. A default scaling of 0.75 is used. Increasing the size of the surfaces usually requires the calculation of further ab initio points as the surface interpolation is more stable for surfaces of limited size. Alternative to the SCALE option, which introduces a uniform scaling of all coordinates, individual scaling of the coordinates as provided by the directive SCALNM may be used.
• SKIP3D=value As the number of 3D and 4D surfaces can increase very rapidly, there exists the possibility to neglect unimportant 3D and 4D surfaces by the keywords SKIP3D and SKIP4D. The criterion for the prescreening of the 3D surfaces is based on the 2D terms and likewise for the 4D terms the 3D surfaces are used. The neglect of 3D surfaces automatically leads to the neglect of 4D surfaces, as the latter depend on the previous ones. By default prescreening is switched on, but can be switched off by SKIP3D=0.0 and SKIP4D=0.0.
• SYM=variable Symmetry within electronic structure calculations can be exploited by the keyword SYM=Auto. Usually this leads to significant time savings. By default this symmetry recognition is switched off as certain calculations may cause some trouble (e.g. local correlation methods). Symmetry in electronic structure calculations may not be mistaken by the symmetry of the mode-coupling terms (see keyword MPG). Once SYM=Auto is used, it is advisable to insert an INT card prior to the call of the Hartree-Fock program.
• THRFIT=value The iterative algorithm for generating potential energy surfaces is based on a successive increase of interpolation points. The iterations are terminated once the interpolation of two subsequent iteration steps became stable. The convergence threshold can be changed by the keyword THRFIT=f. There is currently just one control variable for the different 1D, 2D, 3D, and 4D iterations. The 4 thresholds are different but depend on each other. Consequently, changing the default value (THRFIT=4.0d-2) will change all thresholds simultaneously which keeps the calculation balanced.
• TYPE=variable TYPE=QFF calls a macro, which modifies the parameters of the SURF program in order to compute a quartic force field in the most efficient manner. This implies a reduction of the size of the coupling surfaces and a limitation of the maximum number of points for the $n$D-terms. It should be used for VPT2 calculations. TYPE=ZPVE calls a macro, which changes the defaults for several parameters of the SURF, VSCF and VCI programs. It is meant for the quick and efficient calculation of zero point vibrational energies on cost of some accuracy. For example, the expansion of the potential will be truncated after the 2D terms. As a consequence the output of course is reduced to the presentation of the vibrational ground state only. TYPE=FULL (default) performs a standard calculation as needed for VSCF or VCI calculations. Note that, within XSURF calculations, this keyword will be ignored, but Taylor expansions of the potential can be generated by using the VTAYLOR directive.
• USEMRCC=n Once the Mrcc program of M. Kallay or the Gecco program of A. Köhn is used for determining individual grid points, the option USEMRCC=1 needs to be set, which is needed to ensure proper communication between Molpro and Mrcc. Default: USEMRCC=0.
• VAR1D=variable The SURF program reads the energy of electronic structure calculations from the internal Molpro variables, e.g. ENERGY, EMP2, $\dots$. The internal variable is specified by the keyword VAR1D. Within the example shown above, VAR1D=ENERGY would read the CCSD energy, while VAR1D=EMP2 would read the MP2 energy, which is a byproduct of the CCSD calculation. The default for the VAR1D keyword is the internal variable ENERGY.
• VRC=n Once the keyword VRC=1 is provided, the SURF program will also compute the vibrational-rotational coupling surfaces and thus increases the number of degrees of freedom to 3N-3. Vibrational-rotational coupling surfaces can only be used within the PESTRANS program (see below), but will be neglected in any VSCF or VCI calculations.

The following example shows the input of a calculation which computes energy and dipole surfaces at the MP2/cc-pVTZ level and subsequently determines the anharmonic frequencies at the VSCF and VCI levels. Hartree-Fock calculations will not be restarted and the .log-file is directed to the scratch directory as defined by the $TMPDIR variable.

memory,20,m
orient,mass
geometry={
   3
Water
O          0.0675762564        0.0000000000       -1.3259214590
H         -0.4362118830       -0.7612267436       -1.7014971211
H         -0.4362118830        0.7612267436       -1.7014971211
}

mass,iso
basis=vdz
logfile,scratch

hf
mp2
optg
frequencies,symm=auto

label1
int
{hf
 start,atden}
{mp2
 cphf,1}

{surf,start1D=label1,sym=auto
 intensity,dipole=2}
vscf,combi=1
vci,combi=1

VMULT,options

The level of the electronic structure calculations can be changed for the different $i$-body terms in the expansion of the potential. As a consequence, the keywords START2D, START3D, VAR2D and VAR3D exist in full analogy to the keywords START1D and VAR1D in standard calculations (see above). The number always represents the level of the expansion term. Such calculations are termed multi-level calculations. There does not exist a corresponding set of keywords for the 4-body terms. 4-body terms will always use the variables specified for the 3-body terms (this restriction is lifted in the XSURF program.

MULTI=n The keywords START1D, START2D, START3D in combination with the commands VAR1D, VAR2D and VAR3D allow for the calculation of multi-level potential energy surfaces. This would imply in principle that the 1D term of the potential needs to be computed at all three levels and the 2D term at two computational levels. As certain low level results are byproducts of more sophisticated methods (e.g. the HF energy is a byproduct of an MP2 calculation or the MP2 energy is a byproduct of a CCSD(T) calculation) the computational overhead can be avoided by the MULTI option.
MULTI=1: This is the default and most expensive choice. The 1D potential will be computed at all 3 levels of theory. Likewise, the 2D potential will be calculated at 2 levels explicitly. An example would be:

1D: CCSD(T)/cc-pVTZ
2D: MP4(SDQ)/cc-pVTZ
3D: MP2/cc-pVDZ
{SURF,Start1D=label1
 VMULT,Start2D=label2,Start3D=label3,Multi=1}

MULTI=2: All information is provided by the preceding calculations and thus no part of the potential has to be computed twice. Examples:

1D: CCSD(T)/cc-pVTZ
2D: CCSD(T)/cc-pVTZ
3D: MP2/cc-pVTZ
{SURF,Start1D=label1
 VMULT,Start2D=label1,Start3D=label2
 VMULT,Var3D=EMP2,Multi=2}

1D: CCSD(T)/cc-pVTZ
2D: MP2/cc-pVTZ
3D: MP2/cc-pVTZ
{SURF,Start1D=label1
 VMULT,Start2D=label2,Start3D=label2
 VMULT,Var2D=EMP2,Var3D=EMP2,Multi=2}

MULTI=3: The 2D potential provides all information for the 3D part while there is no connection between 1D and 2D. Consequently, the 1D contributions need to be computed twice (at the 1D and 2D levels) while all other terms will be computed just once. Examples:

1D: CCSD(T)/cc-pVTZ
2D: MP4(SDQ)/cc-pVTZ
3D: MP2/cc-pVTZ
{SURF,Start1D=label1
 VMULT,Start2D=label2,Start3D=label3
 VMULT,Var3D=EMP2,Multi=3}

1D: CCSD(T)/cc-pVTZ
2D: MP4(SDQ)/cc-pVTZ
3D: MP4(SDQ)/cc-pVTZ
{SURF,Start1D=label1
 VMULT,Start2D=label2,Start3D=label2,Multi=3}

MULTI=4: The 1D calculation provides all information for the 2D potential but does not so for the 3D part. Hence, the 1D contribution and the 2D contributions need to be computed twice. Examples:

1D: CCSD(T)/cc-pVTZ
2D: CCSD(T)/cc-pVTZ
3D: MP4(SDQ)/cc-pVTZ
{SURF,Start1D=label1
 VMULT,Start2D=label1,Start3D=label2,Multi=4}

1D: CCSD(T)/cc-pVTZ
2D: MP2/cc-pVTZ
3D: MP2/cc-pVDZ
{SURF,Start1D=label1
 VMULT,Start2D=label2,Start3D=label3
 VMULT,Var2D=EMP2,Multi=4}

In 2D and 4D calculations (i.e. NDIM=2,4) the VMULT command can be used as well. In 4D calculations the last level must always be identical to the 3D level. In 2D the meaning of MULTI=1 and MULTI=3 is the same. Likewise, MULTI=2 and MULTI=4 are the same in case of 2D calculations. START2D=label START2D and START3D define labels in the input stream in order to compute the 2D and 3D terms at different levels of electronic structure theory than the 1D terms. The use of the START2D and START3D commands usually requests the use of GOTO commands in the input. VAR2D=variable The keywords VAR2D and VAR3D are defined in full analogy to the VAR1D option. They specify the internal variable (e.g. ENERGY, EMP2, CCSD, …) to be read out for a given grid point.

The following example shows a 1D:CCSD(T)/cc-pVTZ; 2D:MP4(SDQ)/cc-pVTZ and 3D:MP2/cc-pVTZ multi-level calculation. As the MP2 energy is a byproduct of the CCSD(T) and MP4(SDQ) calculations only the 1D grid points will be computed twice (at the CCSD(T) and MP4(SDQ) levels). The 1D and 2D energies will be obtained from the internal variable ENERGY while the 3D energies make use of the EMP2 variable.

memory,50,m
orient,mass
geometry={
6
Ethene
C          0.0000000000        0.0000000000       -0.6685890718
C          0.0000000000        0.0000000000        0.6685890718
H          0.0000000000       -0.9240027061       -1.2338497710
H          0.0000000000        0.9240027061       -1.2338497710
H          0.0000000000        0.9240027061        1.2338497710
H          0.0000000000       -0.9240027061        1.2338497710
}

mass,iso
basis=vtz
logfile,scratch

hf
ccsd(t)
optg
freq,symm=auto

label1
int
{hf
 start,atden}
ccsd(t)
goto,label4

label2
int
{hf
 start,atden}
{mp4
 notripl}
goto,label4

label3
int
{hf
 start,atden}
mp2

label4
{surf,start1D=label1,sym=auto
 vmult,start2D=label2,start3D=label3,Var3D=EMP2,Multi=3}
vscf
vci

INTENSITY,options

The INTENSITY directive of the SURF program provides the option to alter the electronic structure methods for calculating the dipole surfaces. It also allows to define the VARDIPnD[X,Y,Z] variables separately. $n$ describes the dimension of the coupling surface and can be chosen to be 1 - 4.

Dipole surfaces can be computed for all those methods for which analytical gradients are available in Molpro. For all methods except Hartree-Fock this requires the keyword CPHF,1 after the keyword for the electronic structure method. In multi-level schemes for which the variables VAR1D, VAR2D and VAR3D are set individually, the VARDIPnD[X,Y,Z] variables have to be set accordingly. Symmetry is currently only implemented for the 1D, 2D and 3D dipole surfaces. For 4D terms symmetry will automatically switched off at the moment. The determination of dipole surfaces beyond Hartree-Fock quality effectively doubles the computation time for surface calculations.

• DIPOLE=n Allows to switch between the different dipole surface calculations.=0 switches off all dipole calculations. DIPOLE=1 (this is the default) computes the dipole surfaces at the Hartree Fock level of theory, and therefore does not increase the computation time of electronic structure theory. DIPOLE=2 switches on the dipole surfaces at the full level of theory, therefore CPHF,1 is required. This effectively doubles the computation time for surface calculations.
• NDIMDIP=n This denotes the term after which the $n$-body expansion of the dipole surfaces is truncated. The default is set to 3. Note that NDIMDIP has to be lower or equal to NDIM.
• NDIMPOL=n This variable denotes the term after which the $n$-body expansion of the polarizability tensor surfaces is truncated. The default is set to 2. Note that NDIMPOL has to be lower or equal to NDIM and must be smaller than 4.
• POLAR=n By default (POLAR=0) Raman intensities will not be computed. POLAR=1 switches the calculation of polarizability tensor surfaces on. Note that currently only Hartree-Fock and MP2 polarizabilities are supported, which requires the POLARI keyword in the respective programs. Besides that, the frozen core approximation cannot yet be employed within the calculation of MP2 polarizabilities.
• VARDIP1DX=variable Variable which is used for the $x$ direction of the dipole moment for 1D surfaces.
• VARDIP1DY=variable Variable which is used for the $y$ direction of the dipole moment for 1D surfaces.
• VARDIP1DZ=variable Variable which is used for the $z$ direction of the dipole moment for 1D surfaces.
• VARPOL1DXX=variable Variable which is used for the $xx$ component of the polarizability tensor for 1D surfaces.
• VARPOL1DYY=variable Variable which is used for the $yy$ component of the polarizability tensor for 1D surfaces.
• VARPOL1DZZ=variable Variable which is used for the $zz$ component of the polarizability tensor for 1D surfaces.
• VARPOL1DXY=variable Variable which is used for the $xy$ component of the polarizability tensor for 1D surfaces.
• VARPOL1DXZ=variable Variable which is used for the $xz$ component of the polarizability tensor for 1D surfaces.
• VARPOL1DYZ=variable Variable which is used for the $yz$ component of the polarizability tensor for 1D surfaces.

The higher order terms VARDIPnD[X,Y,Z] and VARPOLnD[XX,$\dots$,YZ] can be defined the same way. An example for a calculation, which provides both, infrared and Raman intensities, is given below.

label1
{hf
 start,atden}
{mp2
 core,0
 polari}

{surf,start1D=label1,ndim=3,info=1
 intensity,polar=1,ndimpol=3
 scalnm,auto=on }

poly,polar=1,ndimpol=3
vscf,pot=poly,polar=1,ndimpol=3,info=1
vci,pot=poly,version=4,polar=1,ndimpol=3,info=1

ALTER,options

The ALTER directive of the SURF program allows to apply error correction schemes for individual single point calculations. For example, in case that the Hartree Fock calculation for a certain grid point did not converge and the ORBITAL directive in the subsequent electron correlation calculation uses the IGNORE_ERROR option, an alternative calculation scheme can be provided, e.g. MCSCF in contrast to RHF. In the case of multi level calculations the ALT2D and ALT3D options can be set according to the START2D and START3D options. Note that the energy variable has to be the same in the original method and the alternative. Within the XSURF program the ALTER directive is defined in a different way.

• ALT1D=label Alternative procedure to calculate the single points of the 1st level.
• ALT2D=label Alternative procedure to calculate the single points of the 2nd level.
• ALT3D=label Alternative procedure to calculate the single points of the 3rd level.

Note that DFT calculations often still converge when RHF calculations already fail to do so.

DISK,options

As SURF calculations are very demanding it is highly recommended to dump the grid representation of the potential to disk. This allows for convenient restarts for subsequent vibrational structure calculations. In contrast to previous Molpro releases, the current version supports only restarts from external ASCII-files - due to many new options. Note, that restarts can also be performed for incomplete potentials. An open-access database of high-level potential files is available at http://http://pes-database.theochem.uni-stuttgart.de/surfaces. Potential files have their own version number. Potential file versions older than 3.0 may not support the full capabilities of the current Molpro release.

• WHERE=value In combination with the keywords DUMP and EXTERN for an external restart file, the keyword WHERE specifies the path for the external ASCII file. Two options are available, WHERE=home and WHERE=scr. As the external files can be huge for SURF calculations, they will be stored on the scratch disk given by the Molpro variable $TMPDIR by default.
• DUMP=file name The potential can be dumped into an external ASCII-file which can be used for restarting. Its name must be provided as the argument of the DUMP keyword, e.g. DUMP=’formaldehyde.pot’. The ASCII-file provides the interface to other programs and offers the possibility for controlled storage and modification of the computed potentials. Dipole and polarizability surfaces will also be dumped if available.
• EXTERN=file name In principle, SURF calculations should be restartable at any point of a truncated calculation. However, since only fully converged difference potentials will be stored in the ASCII file (DUMP), this is not the case. As a consequence, SURF calculations can be restarted from dump-files once a batch of surfaces has been dumped. Since the generation of the 2D or 3D surfaces usually requires about 2 to 3 batches, there are about 6-10 restart points for surface calculation including 3D potentials (see also the keywords BATCH3D and BATCH4D of the SURF program). Restarting from the ASCII dump-file is possible for any type of VMULT calculation. Harmonic frequencies should not be recalculated for restarts.
• SAVE=record Once a complete surface has been generated, a record can be specified, where to dump the potential in the temporary binary files. Once this keyword is not provided explicitly, the potential will be dumped in record 5600.2

LINCOMB,options

The LINCOMB directive allows for the calculation of linear combinations of normal coordinates for the expansion of the potential. This is realized by 2×2 Jacobi rotations. At most 3N-6/2 rotations can be provided in the input.

• ANGLE=value Rotation angle in degree.
• LOCAL=n LOCAL=1 localizes the normal coordinates of the CH-stretchings. Note that this destroys symmetry of these modes. Usually localization has strong impact on subsequent VSCF calculations. LOCAL=3 localizes the normal coordinates of a molecular cluster to the contributing entities. This localization scheme localizes within the individual irreps, which usually leads to a very faint localization. Switching symmetry off by MPG=1 in the SURF program leads to a much stronger localization. LOCAL=2 is a combination of LOCAL=1 and LOCAL=3.
• NM1=n, NM2=m Denotes the normal coordinates to be rotated.
• THRLOC=value (=1.0d-6 Default) Threshold within the localization procedure.

SCALNM,options

The SCALE option of the SURF program enables a modification of the extension of all difference potentials by a common factor. In contrast to that the SCALNM directive allows for the scaling with respect to the individual normal coordinates. This is the recommended choice for potentials dominated by quartic rather than quadratic terms. At most 3N-6 individual scale factors and shift parameters can be provided. In particular the AUTO option was found to be very helpful in practical applications.

• AUTO=on / off AUTO=on (default) switches on an automatic scaling procedure of the potential in order to determine meaningful elongations and SHIFT values with respect to all coordinates, i.e. for each normal mode an optimized scaling parameter SFAC and SHIFT parameter will be determined. Usually this results in an increased number of 1D grid points. The AUTO keyword intrinsically depends on the thresholds and parameters, which can be controlled by the keywords THRSHIFT, ITMAX, LEVMAX, DENSMAX, and DENSMIN.
• DENSMAX=value Threshold for the maximum vibrational density on the edges of the potential needed for the automated upscaling of the potentials (see keyword AUTO).
• DENSMIN=value Threshold for the minimum vibrational density on the edges of the potential needed for the automated downscaling of the potentials (see keyword AUTO).
• ITMAX=n Specifies the maximum number of iterations within the automatic scaling of the potentials (see Keyword AUTO).
• LEVMAX=n Maximum number of vibrational states to be included for controlling the automated scaling and shifting procedure. The default is set to 5. This value should support subsequent VCI calculations.
• MODE=n Denotes the normal coordinate to be scaled or shifted.
• SFAC=value Scaling factor for mode MODE. The default is 1.0.
• SHIFT=n Allows to shift the potential with respect to the specified coordinate by n or -n grid points, respectively. Default: SHIFT=0.
• THRSHIFT=value Threshold controlling the automated shifting of potentials as obtained from the state densities on the lhs and rhs of the potentials. The default is given as THRSHIFT=0.05.

DELETE,surface labels

The DELETE directive allows to eliminate individual surfaces within the multi-mode expansion of the potential. Unlike the SKIP3D and SKIP4D keywords, this directive can only be used once a calculation is restarted from a completed potential energy surface calculation. This directive is meant for studying the impact of individual surfaces or to eleminate troublesome surfaces, which failed to converge in the iterative fitting procedure.

REPAR,options

Within the framework of multi-level calculations (see the directive VMULT), 3D and 4D terms can be modeled. The modeling scheme is based on a reparametrization of the semiempirical AM1 method. Consequently, in the input stream the energy variable to be read in must refer to a semiempirical calculation. After the 2D terms the program optimizes the semiempirical parameters in order to represent the 1D and 2D surfaces best.

• ITMAX1D=n The maximum number of iterations in the local optimization of the semiempirical parameters can be controlled by ITMAX1D and ITMAX2D. The defaults are ITMAX1D=100 and ITMAX2D=150.
• RMS1D=value The keywords RMS1D and RMS2D specify the threshold for terminating the 1D and 2D iterations in the local optimization of the semiempirical parameters. The defaults are given by RMS1D=1.d-6 and RMS2D=1.d-6.

The following example shows the input for a surface calculation in which the 3D terms will be modeled.

memory,20,m
orient,mass
geometry={
   3
Water
O          0.0675762564        0.0000000000       -1.3259214590
H         -0.4362118830       -0.7612267436       -1.7014971211
H         -0.4362118830        0.7612267436       -1.7014971211
}

hf
mp2
optg
freq

label1
abinitio
basis=vdz
int
rhf
mp2
goto,label4

label2
semi,am1
int
rhf

label4
{surf,start1D=label1
 vmult,start2D=label1,start3D=label2,multi=4
 repar}
vscf
vci

CHECK,options

The CHECK directive of the SURF program allows for a quality check of a completed surface. This routine simply computes the exact ab initio energies at randomly selected grid points and compares these values with the interpolated ones, which will be used subsequently for the determination of the wavefunction. This program is fully parallelized.

• LEVEL=n Denotes the level to be checked, i.e. 1 corresponds to 1D, etc. Note, levels below $n$ will not be checked automatically.
• POINTS=value Determines the number of grid points in one dimension to be checked. The default is set to 4.

GRIDCOMP,options

The GRIDCOMP directive of the SURF program allows to interface Molpro with a grid computing client such as Segl. It is also possible to use the grid computing interface without any grid computing client by using the two scripts (create_submit and collect) being supplied in the directory “src/vscf/”. The charge of the molecule as well as some other general commands are transferred to the individual grid point command files, which are printed out in the subdirectory points. If there are any doubts whether the specified command is transferred to the single point or not, the command should be given within the definition of the electronic structure calculations.

• MEMORY=n Denotes the amount of memory (in MW) which is needed in each single point calculation. The default is given by 100 MW.
• MAXFILE=n Defines the maximum number of files produced in a single run of the grid computing interface. Default: MAXFILE=100000.
• FREEZE=n Determines to which permanent file some relevant orbital information has been saved. This is necessary when using the i explicitly correlated or local methods. Default: FREEZE=0.
• FRNAME=value The name of the permanent file, which should be used in each single point.

memory,50,m
orient,mass
geometry={
6
Ethene
C          0.0000000000        0.0000000000       -0.6685890718
C          0.0000000000        0.0000000000        0.6685890718
H          0.0000000000       -0.9240027061       -1.2338497710
H          0.0000000000        0.9240027061       -1.2338497710
H          0.0000000000        0.9240027061        1.2338497710
H          0.0000000000       -0.9240027061        1.2338497710
}

mass,iso
basis=vtz
logfile,scratch

hf
ccsd(t)
optg
freq,symm=auto

label1
int
{hf
 start,atden}
ccsd(t)
goto,label4

label2
int
{hf
 start,atden}
{mp4
 notripl}
goto,label4

label3
int
{hf
 start,atden}
mp2

label4
{surf,start1D=label1,sym=auto
 gridcomp,memory=10
 vmult,start2D=label2,start3D=label3,Var3D=EMP2,Multi=3}
vscf
vci

To generate a potential energy surface with the grid computing interface, follow these steps:

• Run a Molpro calculation on the master control file (see the example above) to generate the control files for the individual single points.
• Calculate the energies for all generated control files in the points subdirectory
• Collect the results and start with point 1 until no more new .com files are produced.

The results of the point calculations should be collected as listed below:

grep -h '*** 1D Surf ilev=1' *.out >> results1Dilev-1
grep -h '*** 1D Surf ilev=2' *.out >> results1Dilev-2
grep -h '*** 1D Surf ilev=3' *.out >> results1Dilev-3
grep -h '*** 2D Surf ilev=2' *.out >> results2Dilev-2
grep -h '*** 2D Surf ilev=3' *.out >> results2Dilev-3
grep -h '*** 3D Surf' *.out >>results3D
grep -h '*** 4D Surf' *.out >>results4D

To reduce the calculation time of the first step, the restart procedure can be used.

It is recommended to

• use the ORIENT,MASS keyword in order to rotate the molecule into standard orientation. This is necessary for a full exploitation of symmetry within the generation of the potential energy surface.
• use the MASS,ISO keyword to use the most available isotopes.
• use multi-level schemes in combination with symmetry and a parallelized Molpro version (MPPX) in order to speed up the calculations. Explicitly correlated methods are preferable over conventional approaches.

Problem: The Surf calculation crashes with an error message like

 ?ERROR IN VIRTORB: INCORRECT NUMBER OF ORB...
 ERROR EXIT
 CURRENT STACK:      MAIN

Solution: The program has problems in the symmetry conversion when restarting a Hartree-Fock calculation from the reference calculation at the equilibrium geometry. You need to start the Hartree-Fock calculations independently by using the keywords start,atden.

Problem: In parallel calculations (mppx) the CPU-time of a SURF calculation differs considerably from the real-time (wallclock time).

Solution: There may be two reasons for this: (1) Usually a SURF calculation spends a significant amount of the total time in the Hartree-Fock program and the 2-electron integrals program. As the integrals are stored on disk, 2 processes on the same machine may write on disk at the same time and thus the calculation time depends to some extend on the disk controller. It is more efficient to stripe several disks and to use several controllers. This problem can be circumvented by distributing the job over several machines, but limiting the number of processors for each machine to 1. (2) The integrals program buffers the integrals. Parallel jobs may require too much memory (factor of 2 plus the shared memory) and thus the integrals buffering will be inefficient. Try to reduce the memory as much as you can. It might be advantageous to separate the memory demanding VCI calculation from the SURF calculation.

XSURF,options [xsurf]

The XSURF program is not just an extension to the old SURF program, but a new program, which works in a completely different manner. However, the syntax for controlling it is very much the same as for the SURF program. In contrast to the SURF program, XSURF can handle $n$-mode and Taylor expansions of the PES of arbitrary order. Moreover, it can handle any kind of symmetry, e.g. non-Abelian point groups or permutational symmetry, Besides that it offers a much more flexible multi-level input and many more options of minor importance. Most importantly, XSURF calculations can be restarted at any point and the external restart files are much smaller and have a completely different structure. In the following, only those options and directives will be listed, which are not valid for the SURF program or which differ considerably. The XSURF program does no longer support a number of keywords being relevant for the SURF program, e.g. BATCH and VRC. The SADDLE keyword is no longer needed as XSURF will recognize automatically if the expansion point is a minimum or a transition state. Moreover, the keywords NDIM, NDIMDIP and NDIMPOL are no longer restricted to values between 1 and 4.
B. Ziegler, G. Rauhut, Rigorous use of symmetry within the construction of multidimensional potential energy surfaces, J. Chem. Phys. 149, 164110 (2018).
B. Ziegler, G. Rauhut, Localized Normal Coordinates in Accurate Vibrational Structure Calculations: Benchmarks for Small Molecules, J. Chem. Theory Comput. 15, 4187 (2019)

• AUTOFIT=n (=0 Default) If AUTOFIT=1, the number of basis function for fitting the grid points is determined automatically. To do so, the fine grid of the energy is compared to the coarse grid points. If the deviation is too high, another basis function is added. The procedure starts with 8 basis functions and stops at the latest at FITXD_MAX basis functions. Once 'AUTOFIT' is used, the keywords FITXD does not have any use.
• CORRECT=n (=1 (on) Default) If a certain subsurface does not converge despite increasing the number of ab initio calculations, symmetry in this subsurface (if any) will be neglected in order to avoid any errors due to inaccuracies in the displacement vectors and the subsurface will be recalculated accordingly. This option is automatically switched off in any Taylor expansions of the PES.
• FITXD_MAX=n(=10 Default) For the automated procedure with AUTOFIT an upper limit for the number of basis functions can be set with this keyword.
• FITMETHOD=n (=1 Default) Within the iterative build-up of the individual subsurfaces, intermediate fitting will be used. This can be based on true multidimensional Kronecker product fitting (FITMETHOD=1) or on fitting along one-dimensional cuts (FITMETHOD=2).
• INFO=n (=1 Default) INFO=0 suppresses any information about the program parameters and symmetry information. INFO=1 refers to the standard output, while INFO=2 provides additional information about the symmetry recognition.
• LOWERDIM=n (=2 Default) In order to determine energy differences for the subsurfaces, contributions of subsurfaces of lower order need to be substracted from those of current order. The lower order energy contribution can be determined in three different ways: LOWERDIM=1 determines these contributions from explicit electronic structure calculations, LOWERDIM=2 uses Kronecker product fitting of the lower order terms and LOWERDIM=3 retrieves this information from the intermediate fine grid representation of the underlying subsurfes.
• MINnD_CRIT1=n Minimum number of ab initio calculations per dimension, which are needed before a subsurfaces can be declared to be converged according to criterion 1. For dimension 1-3 the default is 6, while it is reduced to 4 for all higher order subsurfaces.
• NSA=n (=1 Default) This option prints out some information about the progress of the XSURF calculation. NSA=1 prints this information is an additional file, which will be deleted once the XSURF calculation is completed, NSA=2 prints this information to the console.
• ONLYFREQ=n (=0 (off) Default) If set to 1, the XSURF calculation will be terminated after writing the header of the external restart file, i.e. prior to the calculation of the 1D terms.
• POINT_SCHEME=variable (=NOSHIFT Default) The distribution of ab initio points along a coordinate is determined by a fixed point scheme. This distribution has been generated for potentials, which have not been shifted. For strongly shifted potentials, improved point schemes can be used by the option POINT_SCHEME=SHIFT.
• RDM=n (=0 (off) Default) Degenerate modes can be rotated in a manner, that the corresponding 1D potentials will be identical. By default, this feature is switched off, but can be activated by RDM=1. Typically this results in rotational angles of 45 or 135 degrees.
• RDM_THR=n (=1.0d-10 Default) This threshold controls, if the potentials of two degenerate modes are identical or not. See the keyword RDM=n.
• SKIP=n (=1 (on) Default) By default, (pre)screening of any terms higher than 2D of the PES expansion is switched on. It can be deactivated by SKIP=0.
• SKIPCRIT=n (=1 Default) SKIPCRIT defines the method and thus the criterion for (pre)screening of the high-order terms of the PES expansion. SKIPCRIT=1 activates prescreening and SKIPCRIT=2 screening. In the latter case a label must be defined (see SKIPLABEL).
• SKIPLAB=variable The name of the label in the input stream must be defined, which determines the electronic structure level to be used for screening the high-order terms of a PES expansion.
• SKIP_THRx(i)=value ($x$=1,2 $i$=1,2,…) Definition of thresholds used for prescreening (SKIP_THR1) and screening SKIP_THR2 (see option SKIPCRIT). $i$ specifies the individual orders of the $n$-mode expansion of the PES. The defaults are SKIP_THR1(1)=0.0d0, SKIP_THR1(2)=0.0d0, SKIP_THR1(3)=1.0d-18, SKIP_THR1(4)=1.0d-27 and SKIP_THR2(1)=0.0d0, SKIP_THR2(2)=0.0d0, SKIP_THR2(3)=1.0d-2, SKIP_THR2(4)=1.0d-1. All other thresholds are not yet defined and need to be entered explicitly.
• STOP=i (=0 Default) The calculation of an $n$-mode expansion of a PES can be terminated at lower orders, specified by STOP, while all information for the high-order terms is already provided in the header of the external restart file and the potential information once a multi-level scheme has been used.
• THRDEGx=value ($x$=1,2,…) Threshold used for recognizing degeneracies within the individual orders of the PES expansion. The defaults are THRDEG1=5.0d-5, THRDEG2=1.0d-5,THRDEG3=5.0d-6,THRDEG4=1.0d-6.
• THRFITx=value ($x$=1,2,…) Threshold used for testing on convergence of a subsurface according to criterion 1 for the individual orders of the PES expansion. The defaults are THRFIT1=2.0d-3, THRFIT2=5.0d-3,THRFIT3=5.0d-3,THRFIT4=2.0d-2.
• THRLOWx=value ($x$=1,2,…) Threshold used for testing on convergence of a subsurface according to criterion 2 for the individual orders of the PES expansion. The defaults are THRLOW1=0.0d0, THRLOW2=5.0d-7,THRLOW3=1.0d-6,THRLOW4=2.0d-6.
• THRSED=value (=1.0d-6 Default) Threshold for determining symmetry elements of the molecule.
• THRSYMx=value ($x$=1,2,…) Threshold used for recognizing symmetry within a subsurface of the PES expansion - in dependence on the order of the expansion term. The defaults are THRSYM1=5.0d-5, THRSYM2=1.0d-5,THRSYM3=5.0d-6,THRSYM4=5.0d-6,THRSYM5=1.0d-7.
• DELAUTO=variable(=off Default) If DELAUTO=on, all not converged surfaces of the highest considered dimension are deleted. It only works after a restart from an external potfile.

VIBMODE,options

The VIBMODE directive allows to span the PES only with predefined modes. The following options can be combined in various ways.

• ENERGHIGH=x Modes with a frequency lower than x are used to span the surface (according to the harmonic frequencies)
• ENERGLOW=x Modes with a frequency higher than x are used to span the surface (according to the harmonic frequencies)
• HIGH=n the highest n modes are used to span the surface
• LOW=n the lowest n modes are used to span the surface
• MODE=n Mode which is used to span the surface (can be used multiple times)

GRAPH,options

The GRAPH directive is the interface to programs for visualising potential terms and to provide potential information for any other programs in a most simple manner.

• DIRECTORY=path (=./plots1 Default) Path or name of the directory to be specified for dumping the files for visualisation. See the keyword NDIM.
• GEOM=file name This option dumps the xyz coordinates of all structures used within the electronic structure calculations into an ASCII file.
• MOLDEN=file name This allows to dump a file, which can directly read in by the Molden or Wxmacmolplot programs. This allows for the visualisation of the geometry and the harmonic frequencies of the molecule.
• NDIM=n (=0 Default) This keyword writes all nD surfaces and a corresponding Gnuplot script in a separate subdirectory (plots1) in the home-directory in order to allow for visualization of the computed nD surfaces. E.g. the command “gnuplot plotV1D.gnu” in the plots1 directory produces .eps files for all 1D surfaces. This keyword corresponds to the PLOT keyword of the SURF program.

ALTER,options

As the ALTER directive of the SURF program was slightly confusing, it has been completely redefined for the XSURF program. It allows to apply error correction schemes for individual single point calculations. For example, in case that the Hartree Fock calculation for a certain grid point did not converge and the ORBITAL directive in the subsequent electron correlation calculation uses the IGNORE_ERROR option, an alternative calculation scheme can be provided, e.g. MCSCF in contrast to RHF. The ALTER directive always requests to specify a new label, which replaces the old one. If more than one label shall be replaced, the ALTER directive needs to be called repeatedly.

• NEW=label Specification of the new label.
• OLD=label Specification of the old label.

The following example demonstrates the use of the ALTER directive within a multi-level PES calculation. In case that the RHF calculation for the 1D or 2D terms does not converge, this will be recognized by the UCCSD(T)-F12a calculations. The option orbital,ignore_error=2 prevents a termination of Molpro and the programs tries to enforce convergence by the recipe specified below label5.

memory,500,m
gthresh,optgrad=1.d-7,twoint=1.d-14,prefac=1.d-16
geometry={
O   ,,     -0.4187902305    ,   -2.0024156961    ,    0.0000000000
N   ,,      0.9293931761    ,   -0.1403574221    ,    0.0000000000
N   ,,     -0.3385150234    ,    2.0246339597    ,    0.0000000000
H   ,,     -2.2504006895    ,    1.9819452952    ,    0.0000000000
H   ,,      0.6869538492    ,    3.6185403577    ,    0.0000000000
}

logfile,scratch
basis=vtz-f12
mass,iso

charge=1
{rhf
 start,atden}
uccsd(t)-f12a
optg
freq,symm=auto

{rhf
 start,atden}
{uccsd(t)-f12a,freeze_save=1891.2}

label1
symmetry,auto
charge=1
{rhf
 start,atden}
{uccsd(t)-f12a,freeze_start=1891.2;orbital,ignore_error=2}
goto,label4

label2
symmetry,nosym
charge=1
{uks,b3lyp
 start,atden}
goto,label4

label5
symmetry,nosym
basis=vtz-f12
int
{rhf
 start,atden
 save,2111.2}
{multi
 occ,13
 closed,11
 wf,23,1,1
 canonical,3111.2
 start,2111.2}
{rhf
 start,3111.2}
{uccsd(t)-f12a,freeze_start=1891.2}

label4
{xsurf,sym=auto
 alter,old=label1,new=label5
 vmult,start2D=label1,start3D=label2,multi=4
 disk,where=home,dump='N2H2O.pot'}

VMULT,options

The VMULT directive of the XSURF program offers exactly the same options as within the SURF program. However, if you need a more specialized multi-level scheme, XSURF provides an expert mode. Within this you can for example change the method for 3D and 4D terms. When calculating an ab initio point for an $n$ dimensional body term, two aspects need to be considered. First, you may check whether information for higher-order terms can be retrieved from low-order terms, e.g. the MP2 energy is a byproduct of a CCSD calculation. Second, if the information cannot be retrieved, single-point calculations at different electronic structure levels may be repeated for the low-order terms. All this information is internally stored in a matrix. $$\begin{aligned} M&=&\left(\begin{array}{cccc} 1&0&-1&1\\ 0&1&-2&1\\ 0&0&1&1\\ 0&0&0&1 \end{array}\right)\end{aligned}$$ The columns of the matrix belong to the method and the lines to the dimension. The only important numbers are therefore the upper diagonals. A $1$ means calculate something and a number smaller than one, that the information is somewhere available. This matrix is an example for the following set of methods:

1D: CCSD(T)/cc-pVTZ
2D: CCSD(T)/cc-pVTZ
3D: MP2/cc-pVTZ
4D: HF/minao

For the other multi cases (see above) the matrix elements look like:

• MULTI=1: $$M=\left(\begin{array}{ccc} 1&1&1\\ 0&1&1\\ 0&0&1 \end{array}\right)$$
• MULTI=2: $$M_1=\left(\begin{array}{ccc} 1&0&-1\\ 0&1&-2\\ 0&0&1 \end{array}\right); \qquad M_2=\left(\begin{array}{ccc} 1&-1&0\\ 0&1&0\\ 0&0&1 \end{array}\right)$$
• MULTI=3: $$M_1=\left(\begin{array}{ccc} 1&1&-2\\ 0&1&-2\\ 0&0&1 \end{array}\right); \qquad M_2=\left(\begin{array}{ccc} 1&1&0\\ 0&1&0\\ 0&0&1 \end{array}\right)$$
• MULTI=4: $$M_1=\left(\begin{array}{ccc} 1&0&1\\ 0&1&1\\ 0&0&1 \end{array}\right); \qquad M_2=\left(\begin{array}{ccc} 1&-1&1\\ 0&1&1\\ 0&0&1 \end{array}\right)$$

These expert multi level matrices can be set by MMAT(x,y)=n. Only non-zero elements need to be set.

VADD,options

In practice one may come to the conclusion that the $n$-mode expansion of an existing potential has been truncated too early and one needs to extend the order. Once the level of the electronic structure calculations shall not be changed between the highest existing order and the one to be added, one can simply restart the calculation from the exisiting restart file by increasing the option NDIM in the input stream. However, if one needs to alter the electronic structure level, the VADD directive needs to be used.

• START=label Specification of the label for the order to be added.
• QCOORD=variable (=POT Default) This option specifies from where to take the displacement vectors for the normal coordinates. By default they are taken from the external restart file, but QCOORD=FREQ allows to take them from a new harmonic frequency calculation.

VFREQ,options

Usually, the diplacements vectors of the normal coordinates are retrieved from a preceding harmonic frequency calculation called by the FREQ program. Alternatively, these vectors can be obtained from the XSURF program and the VFREQ directive. However, this alternative is solely based on a twofold numerical differentiation and does not take advantage out of analytical derivatives. However it offers a couple of options, which are not available in the FREQ program.

• COORD=n (=2 Default) Symmetry adapted coordinates, COORD=1, or Cartesian coordinates, COORD=2, may be used within the numerical differentiation.
• METHOD=n (=1 Default) This option specifies the number of points within the numerical differentiation, i.e. METHOD=1 refers to the standard 3-point formula (central differences), METHOD=2 denotes the more accurate 5-point formula and METHOD=3 the 7-point formula.
• PRINT=n (=1 Default) Printout control.
• START=label This sets the label in the input stream to determine the electronic structure level to be used.
• STEP=value (=1.0d-2 Default) This option specifies the step width within the numerical differentiation.

LINCOMB,options

In addition to the directive with the same name in the SURF program, in combination with the XSURF program, the possibility is offered to generate blocks of coordinates, which are then localized independently.

• BLx=n This option assigns mode $n$ to block $x$. For example, the line LINCOMB,BL1=1,BL1=2,BL2=5,BL2=6 for a 4-atomic molecule assigns modes 1 and 2 to a 1st block and thus these two modes will be localized afterwards. Modes 3 and 4 will not be affected and thus refer to standard normal coordinates, while modes 5 and 6 constitute a 2nd block.

VTAYLOR,options

By default, the XSURF program generates an $n$-mode expansion. However, the program structure allows also to retrieve a Taylor expansion of the potential, which is identical with a Taylor expansion obtained by differentiation. In principal Taylor expansions of arbitrary order can be generated, but of course it must be guaranteed that the order of coupling terms is sufficiently high, e.g. a sextic force field cannot be obtained from a NDIM=4 calculation, because this calculation generates coupling terms with at most 4 different indices. In such a case, the missing terms will simply be neglected.

• ORDER=n (=5 Default) Number of basis functions within the Taylor expansion.
• POINTS=n (=5 Default) Number of ab initio points controlling the accuracy of the derivatives (e.g. 5-point formula).
• SCALE=value (=7.0d-2 Default) This keyword controls the step width used and corresponds to the SCALE keyword in the SURF and XSURF programs.
• TYPE=variable TYPE=QFF (corresponds to POINTS=5 and ORDER=5) specifies a full quartic force field.
  

TYPE=SQFF specifies a semi-quartic force field (as used in VPT2 calculations).

TYPE=SEXTIC (corresponds to POINTS=7 and ORDER=7) is the shortcut for a sextic force field.

EXTRADATA,options

The XSURF programs allows to compute energy surfaces, dipole surfaces and polarizability surfaces. In addition to that, arbitrary property surfaces can be generated and dumped into an external restart file.

• NDIM=n Dimension of the $n$-mode expansion to be used for the new property.
• NEL=n Number of data to be read in for one point.
• VARx=variable (x=number) Name of the variable, which shall be read from the input file.

INTERFACE,options

The INTERFACE directive allows for the communication with other programs. It writes information about the individual grid points of a PES to an external ASCII file, which can be processed by other software. Likewise, files in the same structure with additional information from external programs can be read in. After reading in all data points, the date points will be transformed into a fine grid and fitted to polynomials.

• COPY=n (=0 Default) Once new data have been generated in the external ASCII file, the coefficients of the corresponding polynomials can be displayed in the POLY program using the option COEF_INTERFACEx with (x=1,2…). It is also possible to replace the energy or dipole surfaces generated by Molpro by these new quantities by COPY=ENE or COPY=DIP.
• DATA=n (=1 Default) DATA=1 provides detailed information about each single point of the PES in a formatted output. DATA=2 provides the geometry and energy of a given point in a single line. New information about this point needs to be added at the end of the line. DATA=2 prints the displacements along the coordinates and the energy in a single line. Again, new information needs to be added at the end of this line.
• NDIM=n (=0 Default) Dimension of the $n$-mode expansion to which the geometry information shall be dumped.
• NRES=n (=1 Default) Number of columns being added to the external file by an external program.
• SURFACE=n (=0 Default) Information about the energy values printed in the external file. SURFACE=0 refers to absolute energies (minus the reference energy), while SURFACE=1 refers to energy differences belonging to the individual increments of the subsurfaces.
• TYPE=variable (=OUT Default) This option controls, if the file shall be written TYPE=OUT or read in TYPE=IN.
• ZERO=n (=1 Default) If set to 1, geometries of lower orders of the $n$-mode representation will be printed, i.e. the external file contains redundand data. ZERO=0 neglects all redundancies and prints only unique points. As a consequence, an external file generated this way cannot be read in again for technical reasons.
• WFU=file name Specifies the name of the external file.

XGRIDCOMP,options

The grid computing interface of the XSURF program differs considerably from the GRIDCOMP directive of the SURF program and the options do not match. However, the philosophy behind is the same and thus this directive creates batches of single point calculations, which can be processed independently by a Make file generated by the XSURF program.

• CORES=n Number of cores to be used within the grid computing.
• LOOP=n When processing the batches of single point calculations, the grid coumputing interface needs to know, if it is the very first batch (LOOP=0), in which just input files for Molpro will be generated or if it is one of the subsequent batches (LOOP=1), in which also the outputs of processed single point calculations need to be read in.
• MEMORY=n Memory request of the individual single point calculations in MW.
• WFU=file name If additional information need to be read in from a .wfu file, this can be specified here.