From:	Giovanni Scalmani <giovanni "at@at" sg2.csrsrc.mi.cnr.it>
Date:	Tue, 21 Jan 1997 02:21:58 -0800 (PST)
Subject:	Summary: Electrostatic effects in molecular crystals

Dear friends,

some days ago I submitted a question about "Electrostatic effects
in molecular crystals":

>
>I would like to know if someone has already faced
>and (possibly) solved the following problem:
>
>I would like to model - within an ab-initio calculation -
>the crystalline enviroment effect on a single molecule
>in a molecular crystal, i.e. I would like to perform a
>(possibly high level) calculation of a single molecule's
>electronic structure under the influence of a
>neighbouring molecules' charge distribution.
>An iterative procedure should be necessary to attain
>self-consistency between the ab-initio charge distribution
>of the reference molecule and the neighbouring molecules'
>charge distribution, which will be in turn derived from the
>former.
>
>Since I have no deep knowledge of crystallography, I will
>appreciate advises about two issues:
>1) Widely speaking: Is such an approach correct? Can you
>   give me some references in the literature? If I have more
>   than one molecule (or even molecular species) in the cell
>   a similar approach on more than one molecule is feasible,
>   but is it worth while? How should I model the neighbouring
>   molecules' charge distribution (Mulliken charges, charges
>   fitted to the molecular electrostatic potential,
>   Atom-In-Molecule derived charges, orientable atomic
>   dipoles, ... ) ?
>2) More technically: do you know a program which reads the
>   crystalline structure (for example from the Cambridge
>   Database) and performs all coordinates manipulation I
>   will need ? In particular to obtain cartesian coordinates
>   of any specified molecule in the cell surrounded by all the
>   neighbouring molecules whose centres of mass are within
>   a certain radius from the reference molecule's centre of
>   mass ?
>

Here I collect the answers I received. Thank you very much to people
who have contributed.

Giovanni.

=Answer 1=====================================================================

From: R29CLOSE ( ( at ) ) ETSU.ETSU-TN.EDU
Organization: East Tennessee State University

---

  Dear Giovanni:
  I have worked on both these problems.  In particular I wrote a program
to generate a unit cell (or more) from crystallographic data.  I sub-
mitted it to CCL last year as DRAWCRYS.FOR.  I can help you with all
this.  In particular, if you send me your xyz coordinates and the space
group, I can do the calculations.  Then with BABEL I can generate an
ALCHEMY file to view the added molecules.
  The first part is more of a problem.  In MOPAC or AMPAC there is a
description of "sparkles"  These are pseudo charges added to a calculation
to simulate the charge of neighbors.
  But first I need to know what you want to calculate.  Often times
the easiest thing to do is simulate H-bands between neighboring molecules
with water molecules in the right place.  ...

=Answer 2=====================================================================

From: Bart Rousseau 
Organization: University of Antwerp

---

Here are some references on solid state ab initio calculations:

Ab-Initio Studies of Crystal Field Effects in Acetylene.
P.Popelier, A.T.H. Lenstra, C. Van Alsenoy en H.J. Geise
Acta Chemica Scandinavica, A42, 539-543 (1988).

An ab-initio study of crystal field effects : Solid and gasphase
geometry of acetamide.
P. Popelier, A.T.H. Lenstra, C. Van Alsenoy and H.J. Geise
Journal of the American Chemical Society, 111, 5658-5660 (1989).

An Ab-Initio Study of Crystal Field Effects.
Part 3 : Solid and Gas Phase Geometry of Formamide, Modeling the
changes in a peptide group due to hydrogen bonds.
P. Popelier, A.T.H. Lenstra, C. Van Alsenoy, H.J. Geise
Structural Chemistry, 2, 3-9 (1991).

Solids Modelled by crystal Fiels ab-initio methods.
Part 4 : Thermal Vibrational Parameters and Lattice Expansion
Coefficient of the Cubic Phase of Acetylene.
K. Verhulst, A.T.H. Lenstra, C. Van Alsenoy, P. Popelier, H.J. Geise.
Acta Crystallographica, submitted.

Solids Modelled by crystal field ab-initio methods.
Part 5 : The phase transitions in biphenyl from a molecular point
of view.
A.T.H. Lenstra, C. Van Alsenoy, K. Verhulst, H.J. Geise
Acta Crystallographica, B50, 96-106 (1994).

Solid State Moddeling
Part VI : 2,3-diketopiperazine. On the Integration
of crystallographic and spectroscopic evidence.
A.T.H. Lenstra, B. Bracke, B. Van Dijk, S. Maes, C. Van Alsenoy,
Herman O. Desseyn, Spiros P. Perlepes.
Acta Crystallographica, submitted.

Ab-initio Studies of Crystal Field Effects.
Part VII : Structure of 2,3-Diketopiperazine Using a 13-Molecule
Cluster, a Calculation involving 1092 Basis Functions.
Anik Peeters, C. Van Alsenoy, A.T.H. Lenstra, H.J. Geise
International Journal of Quantum Chemistry, 46, 73-80 (1993).

Ab-initio studies of crystal field effects.
Part 8 : Structure of formamide oxime using a 15-molecule cluster.
Anik Peeters, C. Van Alsenoy, A.T.H. Lenstra, H.J. Geise
Journal of Molecular Structure (THEOCHEM), 304, 101-107 (1994).

Solids Modelled by Crystal Field Ab-Initio Methods.
Part 9 : Stereoselective Order-disorder in Tri-ortho-thymotide-
3-Buten-2-ol (2/1) Clathrate.
K. Verhulst, A.T.H. Lenstra, C. Van Alsenoy
Acta Crystallographica, B51, 1016-1020 (1995).

Solids Modelled by Crystal Field Ab-initio Methods.
Part 10 : Structure of alpha-glycine, beta-glycine and
gamma-glycine using a 15-molecule cluster.
A. Peeters, C. Van Alsenoy, A.T.H. Lenstra, H.J. Geise
Journal of Chemical Physics, 103, 6608-6616 (1995).

Solids Modelled by crystal Fiels ab-initio methods.
Part 11 : Integration of Chemical Substitution and Packing Schemes
Exploiting Crystallographic and Spectroscopic Evidence Illustrated
via Acetamide and Thioacetamide.
Koen Verhulst, Stefan Maes, Christian Van Alsenoy, Albert T.H. Lenstra
Journal of Molecular Structure (THEOCHEM), submitted.


=Answer 3=====================================================================

From: Gabriele VALERIO 
Organization: Ecole Nationale Superieure de Chimie de Montpellier

---

Caro Giovanni,

I can give you an answer for your second (technical) question. MOLDRAW
(http://www.ch.unito.it/ch/DipIFM/Software/MOLDRAW/moldraw.html)
is
useful to obtain cartesian coordinates or construct z-matrix from
crystalline structural data.

=Answer 4=====================================================================

From: "Gustavo A. Mercier Jr" 

---

Hi!

I can probably shed some light into your problem.

Your problem is analogous to performing a QM computation of
an active site of a protein including the electrostatic effects
of the surrounding protein "environment".

A lot of work went into this area during the 80's when I was
doing my Ph.D. work. Check out names like Orlando Tapia, Peter
vanDuijnen, Hand Dijkman, Arieh Warshel, Michaeil Levitt,Harel
Weinstein,  and obviously myself.

Warshel and Levitt started the whole thing with a paper
in J. Mol. Biol. in the 70's. There is also an early paper by
Umeyama where he simply included the effect of the "environment"
as a set of point charges in the hamiltonian.

Today you are in luck. Many ab initio packages include the ability
to incorporate "solvent" effects and you can use this feature for
your work. Obviously there are several ways of doing this.

- Reaction Field methods (I believe G9x) as originally implemented
by O. Tapia (Although he did semiempirical methods).

- Multipole Representation. Here you have a set of discrete point charges
centered on the atoms, or even higher multipole (some even centered
on lone pairs etc.). G9x, Hondo, and even ADF allow for computations
with point charges. There is a version of Hondo developed by the
the people at the National Bureau of Standards (USA) that incorporates
higher multipoles. (Check out the work of Morris Krauss).

- Multipole Representation plus polarization. This is probably the most
complete treatment. You can check my work and Hans Dijkman's. We never
implemented this in a very user friendly manner. Indeed, I was doing
serial computations with manual steps in between to do so. I then left
the field. Hans and I developed complementary mathematical treatments.
All of it is predicated in the group function method of McWeeney.

I must point out several caveats. As Hans Dijkman (P. vanDuijnen's grad
student at the time) pointed out, you cannot have a direct interaction
between the "environment" and the SCF density at each SCF step. You need
to have an interaction between the converged density of your molecule
and the converged polarized "environment". For example:

Run 1 : Qm computation in the presence of point charges.
effect 1: polarization of the Qm molecule.

Run 2: Use the converge density of 1 to polarize the environment.
For example, with atom centered polarizabilities.

Run 3: Fix the polarized environment and repeat Run 1 using the new
environment, but starting the SCF from the output of Run 1.
effect 3: Second order polarization of Qm molecule.

Run 4: Repeat 2 with the output of Run 3.

Run 5: Repeat 3 with the polarized environment of Run 4.

 ...

The above is a true implementation of the method of group functions
by McWeeney when one of the groups is "classical" and the other
is Qm. Early work in the area was done by H. Weinstein (my Ph.D.
advisor).

This is tedious and unfortunately I never automated it because
In my case it converged extremely quickly. The polarization of
a charged environment is very small.

The reason I describe all of this is because people took the wrong turn
many times. They simply included the environment during the SCF process
and this led to wrong result in thre presence of extended basis sets
(no convergence!), but "worked" with minimal basis sets. The effect
when you include the environment during the SCF process is to
incorporate instantaneous fluctuations in the density of the
Qm and environment groups. This resulted in spurious term like
a the dispersion term that vanDuijnen described. The densities
at each SCF step are not physical observables that can be used
to compute interactions between groups, as McWeeney shows.

Finally, don't forget the factor of (1/2) that occurs in the energy
term because the cost of polarizing the systems is (1/2) the polarization
energy computed from the converged polarized system (i.e. Qm + environment).
This was also a source of trouble then.

 ...


=Answer 5=====================================================================

From: "Donald E. Williams" 

---

Dear Dr. Scalmani:
	Your questions raise many different issues.  In most of my work I have
considered that a given molecule is not polarized by surrounding molecules in
the crystal.  My procedure is to get a good wavefunction for the gas molecule
(I use gaussian-94 for that), calculate the MEP on a geodesic grid around the
molecule.  Then I fit the MEP, using my program PDM97, with net atomic charges,
if possible, or adding additional sites (e.g. lone pairs) if necessary.
	Then to pursue my interest in ab initio crystal structure prediction I
put one or more molecules in an arbitrary cell and minimize the total energy
including exchange repulsion and dispersion attraction.  This mimization is
carried out with my program mpa/mpg.  If everything works the observed crystal
structure (including space group symmetry) is predicted.
	For your specific interest in locating molecules around a reference
molecule in the crystal, mpa/mpg would of course do that.
	A recent reference can be found in Acta Cryst. 1996, A52, 326-328.
-Donald Williams

=Answer 6=====================================================================

From: iguana-0at0-one.net (Ray Crawford)

---

Giovanni,

        Could you please forward the results of this query to me?  We are
currently trying to do something similar in that we are trying to estimate
the energetic differences between a small molecule crystal conformation and
a bound macromolecule/ligand crystal structures.  The programs we are using
to do this are Gaussian and Spartan.  Just a few considerations that you
might want to keep in mind are:

        1.) crystal structure atomic placements are different than force
field/ab initio atom placements.  What this means is that if you try to do a
Single Point calculation on the single molecule crystal structure, you will
get an extremely high energy as compared to a nearby low energy
conformation.  We found that the majority of this energetic difference in
due to the placement of the H's (both bond lengths (which are VERY
influential) and bond angles/torsions).  Other energetic problems result
from the placements of the heavy atoms (although these are not as influential).

        2.) not many programs allow you to get around these problems.  We
are currently using Gaussian because of it's ablility to use very
explicative z-matrices.  Another option are to use Charmm because of its
"restrant" option.

        I hope these two considerations are useful although they may not be
exactly related/relavent to what you are doing...  But I do realize that we
are examining different aspects of similar problems...  Any info you could
pass along would be greatly appreciated.

=Answer 7=====================================================================

From: "Francois GILARDONI - Gp. J.WEBER - Univ. Geneve" 

(The original message was in Italian. I am sorry to Francois for the errors
 I will put in this translation!)
---

I think that there could be many solutions to your problems:

1. Make use of point charges. Gaussian94 can perform that type of calculations.
   I used instead the DFT-code deMon and I had good results (J.Chem.Phys.,
   104 (19), p.7624). The set of charges can be reduced by using Ewald summation
   technique. In that case, the code I used wasn't able to fit the electrostatic
   field with the desired precision.

2. Make use of Lennard-Jones, Morse or Buckingham functions to fit the potential
   energy curve obtained from ab-initio calculations (involving the nearest
   neighbors). According to my experience, Buckingham functions are the most
   suited ... but it will strongly depend upon the particular situation.
   I tried to model Ru(bz)2 and [Ru(cp)2]2+ crystals with this method.

3. The Car-Parinello method and the program MARVIN (pseudo-potential) could also
   be good choice, but I don't know to much about them.

When you are using point charges, don't forget to use a distribution that
reflects the symmetry of the crystal.


----------------------------------------------------------------------
 ^^^ | SCALMANI Giovanni                 giovanni $#at#$ sg2.csrsrc.mi.cnr.it
 o o | Universita' degli Studi di Milano
  |  | Dipartimento di Chimica Fisica ed Elettrochimica
 \_/ | via C.Golgi, 19                          Phone: ++39-2-26603254
     | 20133 Milano (Italy)                     Fax  : ++39-2-70638129
----------------------------------------------------------------------



Similar Messages
02/21/1996:  Summary: CRYSTAL & all
01/07/1997:  Electrostatic_effects_in_molecular_crystals
01/07/1997:  Electrostatic_effects_in_molecular_crystals
08/01/1996:  Re: CCL:M:Heat of formation calculation using MOPAC.
06/28/1995:  Re:POSTED RESPONSES: Quantitative assessment of novel ligands
10/01/1993:  torsion of conjugated systems -- summary
02/23/1999:  SUMMARY: colour of molecules in the solid state (fwd)
02/16/1997:  Re: CCL:crystal packing calculations and programs
01/08/1997:  Re: CCL:Electrostatic_effects_in_molecular_crystals
08/23/1997:  	Re: CCL:MM parameter development


Raw Message Text