[molpro-user] domains in local calculations

[Ricardo Mata]

Dear Ben,

Your observations are within expected. The combination of Pipek-Mezey
(PM) orbitals with the Boughton-Pulay (BP) criteria for the selection
of domains stumbles upon some difficulties in some specific (and
partly documented) cases:

1- when the basis includes diffuse functions (such as aug-cc-pVXZ)
2- in defining domains for aromatic molecules
3- and yes, they can change due to conformational changes

Number 1 and 2, as you pointed out, are usually solved by an adequate
choice of DOMSEL (alias THRBP) and by deleting some functions in the
localization procedure through the use of CPL,DELETE. In the most
recent Molpro version, the default values for DOMSEL are set according
to the basis set in use. This can be somewhat cumbersome, but the
PM/BP procedure is not a black box solution, and requires some fine
tuning before one dives into more computationally intensive tasks.
Problem 3 is a natural consequence of any domain selection procedure,
at some point the bonds do have to break/form, so changes in the
domains due to changes in the geometry are to be expected. This is one
of the most often criticized aspects of local methods based on the
Pulay/Saebo ansatz.


Now, answering your specific questions:

>Is it just a coincidence that the default options are not sufficient with the avtz basis here, or perhaps a coincidence that they are sufficient with
the avqz basis in this case?

That's the glass half full, glass half empty question. Very small
changes tell the difference between the PM/BP way picking the O as
lone domain atom, or adding the H (you can probably see this in the
criteria and the charge) for the OH bond orbital. One should simply be
aware of these pitfalls. Answering that the problem lies on avtz would
be the same as saying that avtz is more problem-prone than avqz. This
is not necessarily true. That depends on the molecule, conformation,
and so on...

>Is this a symptom that I am making some other mistake, or is it "normal" behavior when using the local methods?

This is perfectly normal, I've seen water with avtz going often awry.
Small changes in the basis or defaults can lead back to the "normal"
domains.

There is another procedure you can use for determining domains which
is based on the Natural Bond Orbital method of Weinhold and coworkers.
Instead of PM you can use Natural Localized Molecular Orbitals (NLMO)
and there is a domain criteria based on the NAO population analysis.
You can combine them with PM and BP, for example:

lmp2,loc_method=nbo                           NLMO/NAO
lmp2,loc_method=nbo,thrbp=0.98        NLMO/BP
lmp2,loc_method=pipek                         PM/BP
lmp2,loc_method=pipek,npasel=0.03    PM/NAO

For more information: Mol. Phys. 125, 2753 (2007).

The NLMO/NAO combination seems to be much more stable than PM/BP, in
the case of water, for example, you should always get the same
domains, independent of the basis set, without ever touching the
npasel option (we recommend a value around 0.03). It also works quite
well for a variety of organic molecules. The only problem is that no
analytic gradients are available. If you are interested in optimizing
structures with LMP2, I would recommend using PM/NAO.

If you have any further questions, feel free to write.

Best wishes,
Ricardo


On Thu, Jun 25, 2009 at 10:26 PM, Ben Shepler<bcshepl at emory.edu> wrote:
> Dear Molpro Users,
> I have been running local MP2 and CCSD(T) calculations on water clusters
> with various correlation consistent basis sets (aug-cc-pVnZ and cc-pVnZ with
> n=2,3,4).  I have been seeing what seems to be an unusual result in the
> selection of the orbital domains.
> The orbital domains I would naively expect to get are as follows:  for each
> water molecule 2 domains consisting of only basis functions from the O atom
> (lone pairs) and 2 domains each consisting of basis functions from the O
> atom and one of the Hs (OH bonds).  For almost all choices of basis set and
> using the same default value of THRBP=0.98 this is in fact what I am
> getting.  For example, here are the domains I get for the water trimer using
> the avqz basis and THRBP=0.98:
>  Orbital domains
>    Orb.   Atom     Charge      Crit.
>    4.1    3 O7      1.12       0.000
>           8 H8      0.60       1.000
>    5.1    3 O7      1.10       0.000
>           9 H9      0.56       0.999
>    6.1    3 O7      1.34       0.997
>    7.1    3 O7      1.26       0.999
>    8.1    2 O4      1.26       0.999
>    9.1    2 O4      1.34       0.997
>   10.1    2 O4      1.12       0.000
>           6 H5      0.59       1.000
>   11.1    2 O4      1.10       0.000
>           7 H6      0.55       0.999
>   12.1    1 O1      1.10       0.000
>           5 H3      0.56       0.999
>   13.1    1 O1      1.34       0.997
>   14.1    1 O1      1.26       0.999
>   15.1    1 O1      1.12       0.000
>           4 H2      0.60       1.000
>
> I get the same domains for the vdz, vtz, vqz, avdz, and avqz basis sets.
>  The lone exception is the avtz basis where for the water trimer I get the
> domains:
>  Orbital domains
>    Orb.   Atom     Charge      Crit.
>    4.1    3 O7      1.26       0.000
>           8 H8      0.53       1.000
>    5.1    3 O7      1.54       0.998
>    6.1    3 O7      1.28       0.989
>    7.1    3 O7      1.30       0.999
>    8.1    2 O4      1.30       0.999
>    9.1    2 O4      1.53       0.997
>   10.1    2 O4      1.26       0.000
>           6 H5      0.53       1.000
>   11.1    2 O4      1.28       0.988
>   12.1    1 O1      1.28       0.989
>   13.1    1 O1      1.53       0.998
>   14.1    1 O1      1.30       0.999
>   15.1    1 O1      1.26       0.000
>           4 H2      0.53       1.000
> Further, if the avtz basis is used on the O atom and only vtz is used for
> the hydrogens I get the "expected" domains.  I can also get the "expected"
> domains with avtz by increasing THRBP to 0.99, or by using the option
> pipek,delete=2.
> I see the same behavior in the water dimer, trimer, tetramer and hexamer.
> I understand that it is important to increase THRBP as the basis set is
> increased, and that diffuse functions can cause problems in the localization
> procedure.  Therefore, on its own I am not surprised that the default
> options and the avtz basis do not give the "expected" domains.  It is a
> little surprising in light of the fact that the default options and the avqz
> basis set do.  I was also surprised to see the same behavior for the several
> different clusters.  Furthermore, when using the avtz basis the domains will
> change when looking at different configurations of the clusters.  Again,
> this is perhaps not surprising on its own, but it has not been the case for
> any of the other basis sets I have tested.
> Is it just a coincidence that the default options are not sufficient with
> the avtz basis here, or perhaps a coincidence that they are sufficient with
> the avqz basis in this case?  Is this a symptom that I am making some other
> mistake, or is it "normal" behavior when using the local methods?
> Any insight anyone can share would be greatly appreciated.
> regards,
> Ben Shepler
> _______________________________________________
> Molpro-user mailing list
> Molpro-user at molpro.net
> http://www.molpro.net/mailman/listinfo/molpro-user
>
>