Subject: Re: CCL:M:Questions re dielec. const. & conformers

From: [email-masked] (Peter Shenkin)
Date: Fri, 1 Sep 95 12:29:17 -0400

> From [email-masked]  Fri Sep  1 08:56:02 1995
> 
>      I would like to be able to locate global minimum energy conformations
> for molecules not only in the gas phase but also in solvents with different
> dielectric constants, e.g., n-octanol or water.  I have examined several
> molecular mechanics programmes, which have an adjustable dielectric constant
> parameter as well as an option for using a distance-dependent dielectric.
> However, the literature accompanying the programmes do not discuss at length
> the matter of adjusting the dielectric....

BatchMin, which is the part of the MacroModel package that performs
energy calculations, employs an implicit solvation model that would
appear, in principle, to do what you want.  We utilize the GB/SA
model, which uses an approach based on the generalized Born equation
for electrostatics (see JACS, 1990, 112, 6127-6129) and a conventional
area-based method of accounting for the "solvophobic" effect (generalization
of "hydrophobic").  For an example of its application, see:

  McDonald, D. Q.; Still, W. C., 
  "Conformational Free Energies from Simulation: Stochastic
  Dynamics/Monte Carlo Simulations of a Homologous Series of 
  Gellman's Diamides", 
  J. Am. Chem.  Soc. , 1994, 116,11550.

The problem with using just a dielectric constant for electrostatics 
is that this doesn't account for the displacement of some volume of 
the solvent by the molecule itself.  For example,  if one is
calculating the interaction of, say, a charged pair of atoms from
the ends of glycine zwitterion, one has to account for the fact
that part of the solvent between them is displaced by the intervening
atoms of the GLY (for example, the carbon-alpha).  These displaced
regions can be viewed as regions of low dielectric embedded in a
possibly high-dielectric solvent.  Clearly, the precise geometry
of the situation should ideally be taken into account; electrostatic
theory demonstrates that regions of displaced solvent not "between"
but "outside of" the interacting atoms have to be taken into account
as well.

The effects of the shape of the molecule can be taken into account
reasonably rigorously, given a full charge distribution, by a
solution of the Poisson-Boltzmann equation.  This is carried out
by the well-known program DELPHI and the somewhat less well-known
GB-SOLV.  However, this rigorous computation is too slow to
be used on the fly during molecular mechanics simulations.

The GB/SA model is an attempt to approximate these effects, and
is fast enough to be used in such simulations.  We supply 
parameterizations for water and chloroform right now.

Other methods are used as well.  Purely area-based models have 
been proposed by Eisenberg and by Scheraga;  they take only the 
local environment into account.

A still simpler attempt to incorporate these effects involves
the use of a distance-dependent dielectric;  this doesn't even
take local environment into account.

>      (i) For molecular mechanics calculations, what is the _cut-off_ value
> for the dielectric constant above which the results can no longer be regarded
> as reliable?  5?  10? 25?  Does this upper limit for the value of the di-
> electric constant depend on the force field implemented, or is it applicable
> to molecular mechanics calculations employing _any_ force field?

Perhaps I misunderstand you, but if not, I think you misunderstand the 
meaning of the term "cutoff".  These values (5, 10, 20) are in 
*Angstroms*, not in the units of dielectric constant.  Atom pairs 
further apart than this cutoff are considered not to interact.  
Conventionally used values for this cutoff are in the range, say, 
12 to 20 Angstroms, for electrostatics.

>     (ii).... 
> ...Can anyone offer any practical advice regarding the use of the
> distance-dependent dielectric option in molecular mechanics calculations of
> molecules in organic solvents?  Of course, references to any publications
> directly addressing this matter are most welcome.

Assuming that you really are talking about cutoffs in Angstroms,
we can address this algebraically.  Suppose Rw is the cutoff used
for water;  call it 12 for now.  We want to calculate Ro, the
cutoff to be used for some organic solvent.  Let's assume the
organic solvent has a dielectric constant (Do) of 5 for illustration,
and that water has a dielectric constant (Dw) of 80.  If the
energy of a charged pair at the cutoff is considered too small
to be considerable, then this energy for some charged pair in
water is:

	E_cutoff,w = q1*q2/(Dw*Rw**2)

Now, for the same charged pair in the organic solvent, we want
to pick the cutoff so that E_cutoff has the same (inconsiderable)
value:

	E_cutoff,o = Ecutoff,w = q1*q2/(Do*Ro**2)

This gives:  (Ro/Rw)**2 = Dw/Do, or (Ro/Rw) = SQRT( Dw/Do )
Using the above values for Dw and Do, we get (Ro/Rw) = 4.
So if Rw = 12 in water, you have to use Ro = 36 in an
organic solvent with a dielectric constant of 5.  

A 36-Angstrom cutoff sounds like a huge number, but most
molecules studied in organic solvents are small anyway, so
you might as well run with the complete pair-list.  Incidentally,
for small to medium-sized systems, this is what we recommend
for the GB/SA model.  For GB/SA, long cutoffs are especially
helpful, since the model takes into account the entire 
environment, not just the interacting pair, and not even just
this pair plus its neighboring atoms.

>    (iii) In relation to the preceding two questions, I have a molecular
> mechanics programme with a "dihedral driver" that will systematically rotate
> a maximum of two bonds in a molecule per run, to generate an m x n matrix
> in which each element corresponds to a particular conformation (and steric
> energy) generated during that run.  Since the molecule in which I am inter-
> ested has more than two conformationally-interesting bonds, I thought to ap-
> ply the dihedral driver stepwise to pairs of bonds in order to locate a glo-
> bal minimum energy conformation (lowest steric energy conformation of the en-
> tire conformational space sampled) _in vacuo_....
> ...However, I am not confident that the COSMO run will lead
> to a global minimum energy conformation for the molecule in the solvent...

If I understand you right, you're trying to minimize with molecular
mechanics in vaccuo, holding certain dihedral angles constant, at a
number of values of dihedral constant, then feed the input conformations
to MOPAC.  It would seem that if any of the above solvent models have
any validity at all, you'd be better off preparing your input conformations
with the molecular-mechanics solvation model turned on.  Then, if
you try to run MOPAC with solvation turned on (which I understand 
you to be saying -- I've never used MOPAC), hopefully your input
conformations would be closer to the true quantum mechanical solvated
minimum.  I'd think this procedure would speed the quantum calculation,
as well as diminish somewhat the chance of getting stuck in a local
minimum.

Incidentally, you have to be careful if the results of the quantum
calculation are to be used later in further molecular mechanics
runs.  The solvation models that are associated with various
force-fields are parameterized for use with various sorts of
charge sets.  For example, our GB/SA model works best when
the charges are closest to the in-vaccuo (not the solvated)
fitted charges from quantum calculations.  This is at least
in part because of the way GB/SA has been parameterized.

This has been wordy, but I hope it helps.

	-P.
************ "There Won't Be Any More."  Charlie Rich, RIP ***********
*Peter S. Shenkin, Box 768 Havemeyer Hall, Chemistry, Columbia Univ.,*
*NY, NY  10027;  [email-masked];  (212)854-5143;  FAX: 678-9039*
********* John Gilmore, RIP. Not to mention S. Chandrasekhar. *********