[MMTK] bulding water ice; problems with Peter's code

Konrad Hinsen hinsen@cnrs-orleans.fr
Mon, 14 Jun 1999 11:22:46 +0200

> Unfortunately, Peter's scripts don't seem to run under the new version 
> anymore. I've made some trivial fixes (reference to Amber94ForceField) 
> but I'm stuck with the following error message:
> eugene.leitl@lrz:~/download/mmtk/MMTK-2.0a5/Examples/peter > python hello.py
> Traceback (innermost last):
>   File "hello.py", line 5, in ?
>     block1 = Group('diamond111')
> NameError: Group

You can import Group from MMTK.ChemicalObjects. However, I don't see
why you would need it to build an ice crystal. The class Group
represents chemical groups for use in building molecules.

> Also, I'm somewhat overwhelmed with Lattice.py. I have the following
> data on water ice crystal cell:
> space group I_h (no. 194 P6_{3/mmc})
> unit cell constants:
> a=b=4.516 Angstrom
> c=  7.354 Angstrom
> alpha=beta= 90 deg
> gamma     = 120 deg
> fractional coordinates:
> O  (0.3333, 0.6667, 0.0629)
> H_a(0.3333, 0.6667, 0.1989)
> H_b(0.4551, 0.9102, 0.0182)
> How do I use BravaisLattice for this?

I am not at all a crystal expert, and consequently crystal support
in MMTK may be unconventional or even insufficient. BravaisLattice
needs the three vectors defining the elementary cell; these can
of course be determined from the crystal parameters you have.

The following code should help you with your crystal construction. It
does in fact produce the right ice crystal, but the construction of
the elementary vectors is "specialized" for beta=90 degrees (I don't
remember the general algorithm, and I am too lazy to derive it...).

from MMTK import *
from MMTK.Lattice import BravaisLattice

a = 4.516*Units.Ang
b = 4.516*Units.Ang
c = 7.354*Units.Ang
alpha = 90*Units.deg
beta = 90*Units.deg
gamma = 120*Units.deg

ex = Vector(1., 0., 0.)
ey = Vector(0., 1., 0.)
ez = Vector(0., 0., 1.)

# Caution: this works only for beta=90*Units.deg
a1 = a*ex
a2 = b*Rotation(ez, alpha)(ex)
a3 = c*Rotation(ey, gamma)(ex)

water = Molecule('water')
water.O.setPosition(0.3333*a1 +  0.6667*a2 + 0.0629*a3)
water.H1.setPosition(0.3333*a1 + 0.6667*a2 + 0.1989*a3)
water.H2.setPosition(0.4551*a1 + 0.9102*a2 + 0.0182*a3)

universe = InfiniteUniverse()
for point in BravaisLattice((a1, a2, a3), (3, 3, 3)):
    molecule = copy(water)

> Also, I am new to the Amber94 force field, having only used CHARMm
> before (X-PLOR, EGO). Are there any dramatic differences?

No. The functional form is almost the same, but the parameters were
fitted using different techniques and different criteria. My
impression is that people choose one or the other for political or
technical reasons rather than for scientific ones. For example, I
picked Amber for MMTK because it is freely available and documented on
the Web.

Occasionally I get requests for implementing CHARMM. This is not at
all difficult in MMTK 2.0, but it requires a precise description of
the CHARMM force field, i.e. a complete recipe for finding all energy
terms for a given system and assigning the correct parameters. To my
knowledge, the only complete description is the CHARMM source code,
which I don't have and don't want to look at.

> What are your experiences with MMTK's numerical performance, btw?

Very machine-dependent. The problem is that many C compilers do not
optimize numerical code well, so Fortran programs often have an
advantage. As an example, I recently did some test runs on an IBM
Power3 processor (200 MHz), and found that it is somewhat slower than
my Pentium II (450 MHz) for MD with MMTK, but faster by a factor 2 for
most Fortran code.

Konrad Hinsen                            | E-Mail: hinsen@cnrs-orleans.fr
Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-
Rue Charles Sadron                       | Fax:  +33-
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais