Properties and expectation values

The property program allows the evaluation of one-electron operators and expectation values. Normally, the operators are computed automatically when using the global GEXPEC directive (see section One-electron operators and expectation values (GEXPEC)) or the EXPEC or TRAN commands in the SCF, MCSCF, and CI programs. The explicit use of the property program is only necessary in the rare case that the user is interested in an orbital analysis of the properties.

PROPERTY

invokes the property program.

DENSITY [,record.file] [,specifications]

If this card is present, the density matrix will be read from record record.file and property expectation values will be calculated. If the specification record.file is omitted, the last dump record is used. Density matrices for specific states can be selected using specifications, as explained in section selecting orbitals and density matrices (ORBITAL, DENSITY). Note that the density matrices are stored in the same record as the orbitals.

ORBITAL [,record.file]  [,specifications]

If this card is present, the orbitals are read from record record.file and an orbital analysis of the expectation values is printed (the density matrix must also be provided!). If record.file is omitted, the last dump record is used. This is only meaningful for diagonal density matrices (SCF or natural orbitals). Natural orbitals for specific states can be selected using specifications, as explained in section selecting orbitals and density matrices (ORBITAL, DENSITY).

The required operators are specified by code words. Optionally, the geometry or the nuclear centre at which the operator is computed can be specified.

For each operator, an input card of the following form is required:

code,centre,x,y,z,,factor

code specifies the property. The available operators are given in section One-electron operators and expectation values (GEXPEC).

The other parameters have the following meaning:

• centre row number of Z–matrix or atomic symbol defining the centre at which property shall be calculated; if centre$\ne 0$ you need not read in coordinates.
• x,y,z cartesian coordinates of the point (only if centre=0).
• factor the operator is multiplied by this factor. The default is factor=1 except for REL. In this cases proper factors for relativistic corrections are used unless factor is given. The two commas before factor are needed to preserve compatibility with Molpro96.

PRINT,print

This card is used to control output, mainly for debugging purposes.

• print$= 0$ no test output (default)
• print$> 0$ operators are printed.

The following example computes the dipole quadrupole moments of water and prints an orbital analysis. By default, the origin is at the centre of mass, and this is taken as origin for the quadrupole moments.

examples/h2o_property.inp

***,h2o properties
geometry={o;h1,o,r;h2,o,r,h1,theta}   !Z-matrix geometry input
r=1 ang                               !bond length
theta=104                             !bond angle
hf                                    !do scf calculation
property                              !call property program
orbital                               !read scf orbitals
density                               !read scf density matrix
dm                                    !compute dipole moments and print orbital contributions
qm                                    !compute quadrupole moments and print orbital contributions
{multi;wf;state,2;dm                  !do full-valence CASSCF
natorb,state=1.1                      !compute natural orbitals for state 1.1
natorb,state=2.1}                     !compute natural orbitals for state 2.1

{property                             !call property program
orbital,state=1.1                     !read casscf natural orbitals for state 1.1
density,state=1.1                     !read casscf density matrix for state 1.1
dm                                    !compute dipole moments and print orbital contributions
qm}                                    !compute quadrupole moments and print orbital contributions

{property                             !call property program
orbital,state=2.1                     !read casscf natural orbitals for state 2.1
density,state=2.1                     !read casscf density matrix for state 2.1
dm                                    !compute dipole moments and print orbital contributions
qm}                                   !compute quadrupole moments and print orbital contributions

Alternatively, the dipole and quadrupole moments can be computed directly in the SCF and MCSCF programs, but in this case no orbital contributions are printed:

examples/h2o_gexpec1.inp

***,h2o properties
geometry={o;h1,o,r;h2,o,r,h1,theta}   !Z-matrix geometry input
r=1 ang                               !bond length
theta=104                             !bond angle
gexpec,dm,qm                          !global request of dipole and quadrupole moments
hf                                    !do scf calculation
{multi;wf;state,2                     !do full-valence CASSCF
natorb,state=1.1                      !compute natural orbitals for state 1.1
natorb,state=2.1}                     !compute natural orbitals for state 2.1

There is a another example described in the Program Control section, as well as:

examples/h2o_prop.inp

***,H2O properties
r=1.85                                            !define initial bond distance
theta=104                                         !define initial bond angle (=2 theta)
geometry={o;h1,o,r;h2,o,r,h1,theta}               !geometry input
Basis={
s,1,v;c,1.6;                                      !van Duijneveld 13s, first 6 contracted
p,1,v;c,1.4;                                      !van Duijneveld 8p, first 4 contracted
d,1,3d;f,1,2f                                     !Dunning optimized 3d, 2f
s,2,v;c,1.4;p,2,3p}                               !8s3p basis for H

gexpec,qm,ef1,ef2,fg,,0,0,0.5                     !Field gradient at (0,0,0.5)

hf                                                !closed-shell SCF

pop                                               !Mulliken population analysis using scf density

{dma                                              !distributed multipole analysis using scf density
limit,,3};                                        !limit to rank 3

multi                                             !full valence CASSCF

{pop;                                             !Mulliken population analysis using mcscf density
individual}                                       !give occupations of individual basis functions

{dma;                                             !distributed multipole analysis using mcscf density
limit,,3}                                         !limit to rank 3

{ci                                               !MRCI with casscf reference
natorb,2350.2
dm,2350.2}                                        !save MRCI density matrix on record 2350, file 2

{pop;density,2350.2;                              !Mulliken population analysis using MRCI density
individual}                                       !give occupations of individual basis functions

{property;
density,2350.2;                                   !use property program to compute expectation values
orbital,2350.2;
qm;ef,1;ef,2;fg,,,,0.5}                           !same operators as before
---

Any density matrix can be analysed using the distributed multipole analysis described by Stone, Chem. Phys. Lett. 83, 233 (1981). The multipole moments arising from the overlap of each pair of primitives are calculated with respect to the overlap centre, and then shifted to the nearest of a number of multipole sites. By default these comprise all atoms specified in the integral input. However the list of multipole sites can be modified by deleting and/or adding sites, and also by restricting the rank of multipole which may be transferred to any given site. The atomic charges are stored in the MOLPRO variable ATCHARGE. The i’th element in ATCHARGE corresponds to the i’th row of the Z-matrix input.

Options may appear in any order, except DENSITY, which must be first if given.

The present version does not allow generally contracted AO basis sets.

DMA;

This command initializes the DMA program.

DENSITY,record.file [,specifications]

The density matrix to be analysed is that found in record record on file file. If omitted, record.file defaults to current orbital record. If specified, DENSITY must appear first in the input. Density matrices for specific states can be selected using specifications, as explained in section selecting orbitals and density matrices (ORBITAL, DENSITY).

GENERAL;

(default) invokes the normal program, which copes with any geometry.

LINEAR

invokes a faster program which can be used when all the atoms are arranged parallel to the $z$-axis and only the $m=0$ components of the multipoles are required.

LIMIT,name,lmax;

lmax is the highest rank of multipole that is to be calculated by the program. Default (and maximum) is 10 for the general program and 20 for the linear one. If name is specified, the limit applies only to multipole site name.

NONUCLEAR

The nuclear contributions to properties are not to be evaluated.

ADD,name,x,y,z,lmax,radius;

Add a new site at ($x$, $y$, $z$) with the name specified. The multipole rank is limited to lmax if a value is specified, otherwise the value of lmax specified by the LIMIT directive is used. No account is taken of symmetry; every site in a symmetry-equivalent set must be specified explicitly. The radius of the site may also be specified (default 1.0).

DELETE,name

Delete all atoms with the name given from consideration as a multipole site. Note that original atoms from the integral program have names $1,2,3,\ldots$ as printed in integral output. DELETE,ALL deletes all atoms and gives the multipoles with respect to the origin only.

RADIUS,name,r;

Assign radius $r$ to all sites with the name given. The program moves multipoles at an overlap centre $P$ to the site $S$ for which the value of $|P-S|/r(S)$ is smallest. In the absence of a RADIUS directive, all sites are given radius 1.

The multipoles produced by this analysis are given in their spherical harmonic definitions. Explicit formulae for translating between the cartesian and spherical harmonic definitions of the multipole moments are given in, Explicit formulae for the electrostatic energy, forces and torques between a pair of molecules of arbitrary symmetry, S. L. Price, A. J. Stone, and M. Alderton, Molec. Phys., 52, 987 (1984).

For examples of the use of DMA analysis see, Price and Stone, Chem. Phys. Lett. 98, 419 (1983); Buckingham and Fowler, J. Chem. Phys. 79, 6426 (1983).

The following input calculates SCF multipole moments for water.

examples/h2o_dma.inp

***,h2o distributed multipole analysis
geometry={o;h1,o,r;h2,o,r,h1,theta}   !Z-matrix geometry input
r=1 ang                               !bond length
theta=104                             !bond angle
basis=6-311g**
hf                                    !do scf calculation
{dma;limit,,4}                        !results for total multipoles are

POP;

Invokes Mulliken analysis program, which analyses any density matrix into its contributions from s,p,d,f… basis functions on each atom. In the case of Hartree-Fock, density functional, and MCSCF calculations, the density matrix is taken from the last dump record, unless overridden with the DENSITY card. For other methods, the density matrix has to be written and then must be referred to with the DENSITY card.

The subcommands may be abbreviated by the first four characters. The atomic charges are stored in the MOLPRO variable ATCHARGE. The i’th element in ATCHARGE corresponds to the i’th row of the Z-matrix input.

DENSITY,record.file [,specifications]

Take density matrix to be analysed from record record on file file. Density matrices for specific states can be selected using specifications, as explained in section selecting orbitals and density matrices (ORBITAL, DENSITY). Note that the density matrices are stored in the same record as the orbitals.

INDIVIDUAL;

examples/h2o_pop.inp

***,h2o population analysis
geometry={o;h1,o,r;h2,o,r,h1,theta}   !Z-matrix geometry input
r=1 ang                               !bond length
theta=104                             !bond angle
basis=6-311g**
hf                                    !do scf calculation
pop;                                  !Mulliken population analysis using mcscf density
individual                            !give occupations of individual basis functions

If specified, the Mulliken populations of each individual basis function are printed.

The following example shows how to compute the population first on the Hartree-Fock level, and then on the CI level.

examples/ch4_ci_pop.inp

geomtyp=xyz;
geometry={
c            -1.058493    -0.384055     0.000000
h            -1.058493     0.695337     0.000000
h            -0.040825    -0.743829     0.000000
h            -1.567327    -0.743829     0.881326
h            -1.567327    -0.743829    -0.881326
}
basis={
default=def2-SVP
}
hf;
pop;
ci;dm,9000.2;
pop;density,9000.2

NBO,[WITH_CORE=core_option],[LEVEL=level],[KEEP_WBI=wbi_option];

The Natural Bond Orbital Analysis of Weinhold and coworkers (J. Chem. Phys. 83, 735 (1985), J. Chem. Phys. 83, 1736 (1985) and J. Chem. Phys. 78 (1983) 4066) can be called by the use of the NBO card. It reads from a density or orbital record, and performs the necessary transformations to Natural Atomic Orbitals (NAO), Natural Bond Orbitals (NBO) and Natural Localized Molecular Orbitals (NLMO). The latter can also be saved to a record and later used in local correlation treatments (cf. Section PAO-based local correlation treatments). By default, the full orbital space is used. The core orbitals can, however, be left out of the procedure if core_option=0.

One can choose to truncate the transformation series (e.g., only compute the NAO orbitals), with help of the LEVEL keyword. If level=1, only the NAO transformation will be carried out. For level=2 the NBO transformation is performed, and for 3 the NLMO (default).

Sometimes, the NBO procedure will not converge due to a bad ordering on the 2-center bond search. The first run is based on the Wiberg bond index, but the algorithm switches to the atom ordering on the subsequent runs. This can be avoided by the use of the option KEEP_WBI. If wbi_option=1, the Wiberg bond index is used in all iterations.

SAVE,record.file;

The NLMO orbitals are saved in the specified record, together with the NPA charges.

Dipole moments, quadrupole moments etc. and the corresponding polarizabilities can be obtained as energy derivatives by the finite difference approximation. This is most easily done with the DIP, QUAD, or FIELD commands, which add a specified amount of the dipole, quadrupole, or any general, operators to the Hamiltonian. An error will result if the added perturbation is not totally symmetric (symmetry 1). Note that the orbitals must be recomputed before performing a correlation calculation.

DIP,dx,dy,dz;
DIP+,dx,dy,dz;

Add a finite combination of the dipole operators $\vec\mu=(\mu_x, \mu_y, \mu_z)$, $H_1=\vec d \cdot \vec \mu$ ( $\vec d=(\textit{dx},\textit{dy},\textit{dz})$) to the Hamiltonian (both the 1-electron operator and the core energy). DIP+ adds to any existing field, otherwise any previous perturbation is removed.

The perturbed hamiltonian represents a physical system in a uniform electric field with electric field strength $\vec F= -\vec d$. Therefore the corresponding energy-derivative form of the dipole moment projection in this direction can be obtained as $$|\vec F|^{-1}\vec F \cdot \vec \mu = |2\vec d|^{-1}(E(\vec d)-E(-\vec d)) + O(|\vec d|^2)= |\vec d|^{-1}(E(\vec d)-E(\vec 0)) + O(|\vec d|).$$ The diagonal polarisability in this direction can similarly be calculated via $$\alpha_{\vec d, \vec d} = |\vec d|^{-2}(E(\vec d)+E(-\vec d)-2E(\vec 0)) + O(|\vec d|^2).$$

QUAD,qxx,qyy,qzz,qxy,qxz,qyz;
QUAD+,qxx,qyy,qzz,qxy,qxz,qyz;

Exactly as the DIP command, but adds a combination of quadrupole operators to the Hamiltonian.

FIELD,oper1,fac1, oper2,fac2, …;
FIELD+,oper1,fac1, oper2,fac2, …;

Adds one-electron operators oper1, oper2, …with the corresponding factors fac1, fac2, …to the one-electron hamiltonian. The available operators are given in section One-electron operators and expectation values (GEXPEC). An error will result if the added perturbation is not totally symmetric (symmetry 1).

FIELD+ adds to any existing field, otherwise any previous field is removed.

Note that FIELD does currently not modify core polarization potentials (CPP). If CPPs are present, only DIP and QUAD should be used.

The first examples shows various possibilities to add perturbations to the one-electron hamiltonian.

examples/field.inp

***,H2O finite fields
R     =    0.96488518 ANG
THETA =  101.90140469
geometry={H1
          O,H1,R;
          H2,O,R,H1,THETA}
{hf;wf,10,1}           !scf without field

f=0.05

dip,,,f                !add dipole (z) field to h0
hf                     !do scf with modified h0

field,dmz,f            !add dipole (z) field to H0
                       !same result as previous example
hf                     !do scf with modified h0

quad,,,f               !add quadrupole (qmzz) field to h0
hf                     !do scf with modified h0

field,qmzz,f           !add quadrupole (qmzz) field to h0;
                       !same result as previous example
hf                     !do scf with modified h0

field,zz,f,xx,-0.5*f,yy,-0.5*f
                       !add general field; same result as quad above
hf                     !do scf with modified h0

field,zz,f           !same as before with separate field commands
field+,xx,-0.5*f
field+,yy,-0.5*f
hf                     !do scf with modified h0

field                  !remove field
hf                     !scf without field

The second example shows how to compute dipole moments and polarizabilities using finite fields.

examples/h2o_field.inp

***,H2O finite field calculations

r=1.85,theta=104                   !set geometry parameters
geometry={O;                       !z-matrix input
          H1,O,r;
          H2,O,r,H1,theta}
basis=avtz                         !define default basis
field=[0,0.005,-0.005]             !define finite field strengths
$method=[hf,mp4,ccsd(t),casscf,mrci]

k=0
do i=1,#field                      !loop over fields
  dip,,,field(i)                   !add finite field to H
  do m=1,#method                   !loop over methods
    k=k+1
    $method(m)                     !calculate energy
    e(k)=energy                    !save energy
  enddo
enddo

k=0
n=#method
do m=1,#method
  k=k+1
  energ(m)=e(k)
  dipmz(m)=(e(k+n)-e(k+2*n))/(field(2)-field(3))   !dipole moment as first energy derivative
  dpolz(m)=(e(k+n)+e(k+2*n)-2*e(k))/((field(2)-field(1))*(field(3)-field(1)))  !polarizability as second der.
enddo

table,method,energ,dipmz,dpolz
title,results for H2O, r=$R, theta=$theta, basis=$basis
---

Relativistic corrections may be calculated within the Cowan-Griffin approach by computing expectation values of the mass-velocity and 1-electron Darwin integrals; these should be generated using the property integral program with keyword REL The expectation values can be computed within the SCF, MCSCF and CI programs in the usual way using the EXPEC command, again with the keyword REL. The mass-velocity and Darwin terms, and their sum are subsequently available through the Molpro variables MASSV, DARWIN and EREL respectively. An example can be found here.