SCIPCM in G94

[keith (- at -) kbw350.chem.yale.edu (Todd Keith)]

 Re:
 > I have been having problems getting geometries (and in some cases SCF's)
 > to converge when using the scrf=scipcm option in G94. The higher
 > the dielectric constant, the more problem there seems to be. Does
 > anyone have any suggestions on how this problem might be cured? I
 > have already started with prior hessians and the molecules are simple,
 > so the problem seems to be the inclusion of solvent effects.
 >
 > Thanks in advance, John
 >
      This is one of a many such complaints I've seen on the CCL,
 so a reply (possibly long-winded) to the list might be helpful ...
      The problem is probably due to the multi-center surface
 integration procedure which, unfortunately, is used by default for
 SCI-PCM in G94.  Often, a far superior procedure is to use the non-
 default single center surface integration procedure:  when used
 appropriately, it is much faster and the surface integrals much more
 accurate than the default procedure.  Geometry optimizations are much
 more stable using the single center method (a conseqence of the more
 accurate surface integrals).  To invoke the single center integration
 procedure, one must enter some additional input on the input line which
 specifies the solvent's dielectric constant.
 Here is the format and a partial description of the general SCI-PCM
 input line for G94:
       Epsilon Cont NPhi NTheta ISurf
 Epsilon = Dielectric constant of solvent (for example:  78.3 for water).
 Cont = Value of isodensity surface in atomic units.  0.0004 au
        is recommended and the default.  Molecular volumes defined
        by 0.0004 isodensity surfaces consistently agree fairly well
        with values inferred from corresponding measured pure liquid
        molar volumes or partial liquid molar volumes, which suggests
        that solute cavities by such isodensity surfaces are physically
        reasonable.
 NPhi = Number of equally spaced phi values for numerical surface
        integration, or total number of points in a Special grid.
        Number is per origin.  Default is Special 302.
        Legitimate Special grids are currently:  110 146 194 302
        434 590 770 and 974.
 NTheta = Number of Gauss-Legendre theta values for numerical surface
          integration, or 1 to use a Special grid, which is the default.
          If using Gauss-Legendre, this number should be about half of
          NPhi.
 ISurf = Whether to do surface integrals numerically using a single center
         (2) - the center of nuclear charge - or piecewise, using a
         multi-center method (3).  3 is, in principle, more generally
         applicable than 2 and is the default, but (2) is much faster
         and more accurate when applicable, which is often.  (2) is
         applicable when, for every point on the cavity surface, the
         projection of the surface normal vector onto the corresponding
         unit vector from the single center is greater than zero.  The
         SCI-PCM solute geometry optimizations require reasonably
         accurate surface integrals to work well.  The accuracy of the
         SCI-PCM surface integrations can be monitored by the calculated
         total flux of the nuclear electric field thru the cavity
         surface, which should equal 4*Pi*Qnuc, where Qnuc is the
         total nuclear charge within the cavity.  Another measure
         of accuracy is the total surface polarization charge,
         Qpol, as compared to the net charge of the cavity, Qcav,
         which is (1/(4*Pi)) times the net flux of the solute
         electric field through the cavity surface.  If things
         are done correctly then:  Qpol=-(1-1/Epsilon)*Qcav.
      I would strongly recommend using the single center integration
 procedure by default.  As implied above, the G94 SCI-PCM implementation
 calculates useful independent measures of the accuracy of the SCI-PCM
 surface integrals.  Unfortunately, the values of these quantities (along
 with other useful info like the cavity volume and surface area)
 are not printed out by default, but can only be seen (for every SCF
 iteration particularly the final one) by requesting #P and IOp(5/33=2),
 which also causes the A LOT of other printing in the .log file (Fock
 matrices, density matrices, etc.).  Specifying #P alone will cause the
 interesting quantities to be printed out for the first SCF iteration,
 which may be adequate to determine whether the surface integrals
 are sufficiently accurate to proceed.
       So, if one uses the single center procedure, the main question
 is how many integration points to use.  Here is a rough! guide of
 molecular size versus total number of points to use:
 1 to 6 atoms -->  434 points (Special grid:  NPhi = 434, NTheta = 1)
 7 to 11 atoms --> 590 points (Special grid:  NPhi = 590, NTheta = 1)
 12 to 15 atoms --> 770 points (Special grid:  NPhi = 770, NTheta = 1)
 16 to 20 atoms --> 974 points (Special grid:  NPhi = 974, NTheta = 1)
  > 20 atoms --> use equally spaced Phi, Gauss-Legendre theta with
                 the total number of points being at least
                 NAtoms*100 (e.g., for 1800 total points use
                 NPhi = 60 and NTheta = 30).
 Of course, the actual shape of the solute molecule might dictate the
 need for more or less points, or possibly, for very large or
 very complex solutes the complete failure of the single center method.
 Here is an example input file for a HF/6-31G* SCI-PCM geometry
 optimization of N,N-dimethyl acetamide in acetonitrile (Epsilon = 35.9):
 #P HF/6-31+G* Opt SCRF=SCIPCM
 N,N dma HF/6-31+G* opt SCI-PCM e=35.9, Cont=0.0004, 770 pts, single cent
 0 1
 Z-Matrix, molecule specification
 35.9 0.0004 770 1 2
 Note that the use of the default multi-center procedure in the above
 case would require many more integration points (like 1500) and would
 thus be relatively slow and the optimization may not converge due to
 relatively inaccurate surface integrals.
 Good Luck.
 T. A. Keith
 keith (- at -) kbw350.chem.yale.edu