Subject: SCIPCM in G94

From: [email-masked] (Todd Keith)
Date: Wed, 17 Jul 1996 23:14:59 -0400

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
[email-masked]

>From [email-masked] Tue Jul 16 02:39 EDT 1996
Date: Tue, 16 Jul 1996 01:39:48 -0500
From: christin@irene (Christine Byrd)
To: [email-masked]
Subject: MOPAC output files
Content-Length: 1136
Status: RO

I have a few questions regarding Mopac .out and .syb
files.

My observation:
I've noticed that .syb files are created after running
a MOPAC, despite failure of a compound to pass one of
the gradient tests (Herbert's or Peter's).  On the other
hand, .syb files are not generated when MOPAC "abandons"
a calculation for a specific compound.

I was hoping if perhaps someone would explain:

	(1)  The criteria according to which MOPAC chooses
	     to submit a compound to either "Peter's
             Test" or "Herbert's Test."

	(2)  Which consequences not passing one of the
             gradient tests has on the resulting MOPAC
             calculation.

I need to run MOPAC calculations on a rather large database
and would like to run them using SPL Macros, for the sake
of time and accuracy.  Should the failure to pass one of
the tests mentioned above have consequences on the MOPAC
calculation and thus on the .syb file, is there a clever 
way to check for compliance with the tests utilizing a Macro?

I would appreciate any advice.

			   Cheers,

			   Christine Byrd			   Christine Byr			   e-mail:  [email-masked]