Subject: G98 / displacement issue with frozen coordinates

From: Gerald Lushington <[email-masked]>
Date: Thu, 30 Nov 2000 10:18:42 -0500

Dear all,

In working with a colleague on large cluster optimizations of enzyme active
site (via G98), I've encountered a problem that seems to be:

  - part philosophical
  - part methodological
  - part bug (?)

Anyway, the issue has to do with doing geometry optimizations with some degrees
of freedom contrained.  In general it's proving necessary to freeze the 
coordinates for at least a couple of atoms in order to keep the cluster in some
sort of shape comparable to the real biosystem.  The tradeoff, of course, comes
in the fact that forces and displacements will never achieve global minima.
Standard protocol seems to be to run the job until these quantities level off
to some reasonably converged values.  When we apply this strategy we find that
force and net displacement eventually level off, but the net displacement 
remains very high.  This brings me to question number 1:

1)  If we stop the geometry optimization at a structure whose force and 
    displacement values have leveled, but whose displacement is still high,
    would any subsequent vibrational analysis be at all physically meaningful? 
    Are there any effective guidelines for this?

My first thought on this was that Gaussian would likely decompose both force
and displacement as a sum over coordinates.  Since the only freezing we do is
way out on the perifery of the cluster, while the structural and vibrational
information of interest is in the center of the system, maybe we can track the
forces and displacements in the central region of interest and feel comfortable
about the results when those central contributions are low.

So, to try to follow that strategy, we've been refering to a table generated as
standard G98 output during the geometry optimization that looks like:

 Variable       Old X    -DE/DX   Delta X   Delta X   Delta X     New X
                                 (Linear)    (Quad)   (Total)
    X1       -2.50100  -0.00465   0.00000   0.00000   2.50100   0.00000
    Y1        8.90656   0.01749   0.00000   0.00000  -8.90656   0.00000
    Z1       -8.10547  -0.01647   0.00000   0.00000   8.10547   0.00000
   X11       -3.32922   0.00421   0.00000   0.00000   2.50100  -0.82822
   Y11       -9.47808  -0.00542   0.00000   0.00000  -8.90656 -18.38465
   Z11       -1.01768  -0.00913   0.00000   0.00000   8.10547   7.08779
   X38        3.49152   0.00310   0.00000   0.00000   2.50100   5.99252
   Y38       14.83066   0.00184   0.00000   0.00000  -8.90656   5.92410
   Z38       -5.24649  -0.00415   0.00000   0.00000   8.10547   2.85898
   X48       -1.96597   0.00535   0.00000   0.00000   2.50100   0.53502
   Y48       -4.25077  -0.01552   0.00000   0.00000  -8.90656 -13.15733
   Z48        9.58721   0.02164   0.00000   0.00000   8.10547  17.69267
   X50        1.54048  -0.00800   0.00000   0.00000   2.50100   4.04148
   Y50       -8.97461   0.00161   0.00000   0.00000  -8.90656 -17.88117
   Z50        5.17857   0.00810   0.00000   0.00000   8.10547  13.28404
    R1        2.84387   0.01822   0.00000   0.00000   0.00000   2.84388
    R2        2.02201   0.01493   0.00000   0.00000   0.00000   2.02201
    R3        2.02201   0.01403   0.00000   0.00000   0.00000   2.02201
    R4        2.02201   0.01490   0.00000   0.00000   0.00000   2.02201
    R5        2.81303  -0.01054   0.00000  -0.01044  -0.00995   2.80307
    R6        2.10327  -0.02056   0.00000  -0.01706  -0.01706   2.08621
    R7        2.10327  -0.01967   0.00000  -0.01633  -0.01633   2.08694
    R8        2.60262  -0.00578   0.00000  -0.00448  -0.00457   2.59805
    R9        2.55751  -0.00018   0.00000  -0.00230  -0.00185   2.55566
    ..           ..                                                ..

NOTE:  in the above, the first bunch of coordinates (X1-R4 inclusive) are all
       frozen.  R5 onwards (many coordinates not shown) are freely varying.

Anyway, in the variable coordinates toward the bottom, it seems that the 
"DE/DX" is a real force, the "Delta X" are real displacements, and the "Old X" 
and "New X" are exactly as they sound -- initial and final values for the 
coordinate.  However:

2) What's going on with these the "Delta X" and "New X" values for the frozen
   coordinates??  They're frozen (and really do not vary during the 
   optimization) so shouldn't "New X" = "Old X" regardless of nonzero forces 
   and predicted displacements?  

3) Furthermore, the values written in "Delta X (Total)" look extremely fishy to
   me?  They're all either "2.50100", "-8.90656", "8.10547" or "0.00000"!  Is
   this a bug, and if so, might it be corrupting the value given for net
   displacement at the end of each geometry step?

Anyway, this may be a complex issue, but any insight you may have would be
greatly appreciated!  

                                        - Gerry

P.S.  I did drop Gaussian Inc. a note on this several weeks ago -- I guess it's
      somewhere in their queue.

----------------------------------------------------------------------
Gerald H. Lushington               Ph: 614-292-6036
Research Specialist                Fax: 614-292-7168
OSC / PET-CCM                      e-mail: [email-masked]
1224 Kinnear Road                  http://www.arl.hpc.mil/PET/cta/ccm
Columbus OH 43212-1163             http://www.asc.hpc.mil/PET/CCM

---

From: "Stefan Fau" <[email-masked]>
Date: Thu, 30 Nov 2000 15:31:20 -0500

Hi,

an interesting practical problem... Here are my 2 cents,
comments invited.

to 2) and 3)
Delta X and New X are the differences and new values
of the geometry descriptors. G98 allows mixed geometry
specifications using internal coordinates together with
cartesian coordinates. Freezing cartesian coordinates is
equivalent to freezing all internal coordinates between this
group of atoms (provided you have only one spatial group
of frozen atoms).

If you have two clusters of frozen atoms linked by a
unconstrained part of the molecule, one way of optimizing
the molecule is to allow movement of the two clusters with
respect to each other. I don't think G98 does it this way, but
you can check by tracking the distances between atoms of
different frozen clusters over the optimization.

Another reason for blockwise movement of the frozen atom
clusters is that sometimes the standard orientation of the
molecule flips from one optimization step to the next. (I think
that the standard orientation is determined by the moments of
inertia. Sometimes two of them are very similar and change
sequence from one geometry to the next.)

I think the list of geometry variables contains both the original
set of geometry descriptors and the redundant internal coordinates.
I am not sure if original descriptors that are NOT redundant internal
coordinates enter in the computation of the displacement used as a
convergence criterion. You may want to browse the CCL archives:
I faintly remember a mail from Gaussian (maybe Doug Fox) on that
topic a year or two ago.

to 1)
Vibrational frequencies are molecular properties. They all depend
on all of the molecule. At the same time, some are quite specific
for functional groups. Therefore the influence of other parts of the
molecule must be small, probably rapidly declining with distance.

The common implementation of harmonic frequencies requires the
forces (=first derivatives of the energy w.r.t. nuclear positions) to
be
(nearly) zero. I don't know how the second derivatives change with
displacements that cause nearly no forces, but the displacements must
be along floppy modes with very small force constants and frequencies.
I suppose that big residual displacements will mainly influence small
frequencies.

In the harmonic approximation you may diagonalize the matrix of
mass weighted second derivatives of the energy with respect to
nuclear coordinates and get the vibrational modes and vibrational
energies. Like with the solution of the HF equation, you can apply
unitary transformations to get a localized description (ideally only
one internal coordinate varying in a mode) that is not diagonal
anymore. The off-diagonal matrix elements describe how the localized
modes influence each other.

More to your problem: I think you can use frequencies that have NO
movement of frozen atoms in the corresponding vibrations. You may
try to block-diagonalize the matrix with a unitary transformation that
constricts vibrations with motion of frozen atoms in one block. Since
the coupling elements with the other block would be zero, the
frequencies
>from that block would be valid.

One caveat:
When you truncate a molecule, you omit interactions between the cut
off
parts as well as interactions between the cut off and the retained
part of
the molecule. If you have to freeze clusters of atoms to force the
truncated
molecule into the original configuration, some of the omitted
interactions
were important. I would like to make sure that none of the cut off
parts had
a significant interaction with the interesting part of the molecule.

I hope it answers some of your questions,
        Stefan
______________________________________________________________________
Dr. Stefan Fau                    |      [email-masked]
Quantum Theory Project     |     (352) 392-6714
University of Florida
Gainesville, FL 32611-8435