Subject: Re: CCL:G98 / displacement issue with frozen coordinates

From: Gerald Lushington <[email-masked]>
Date: Thu, 7 Dec 2000 14:58:08 -0500

Greetings once again.

I also received a helpful note from Frank Jensen in response to my original
query.  He also sent a program which level-shifts away the effects of the 
frozen atoms which I think is an effective means for addressing my point 1).

Anyway, his note is appended below.  Thanks again to all who responded!

                                   - Gerry

Forwarded note:

Date: Fri, 1 Dec 2000 07:49:46 +0100 (MET)
From: Frank Jensen <[email-masked]>
To: Gerald Lushington <[email-masked]>
Subject: Re: CCL:G98 / displacement issue with frozen coordinates

        Gerald,
        I have encountered the same problem some time ago,
it is clearly a bug in G98. I reported it to Gaussian,
and they put it on the 'fix'-list. This was with rev. A7,
I don't know if it has been fixed in rev. A9. Anyway,
once the forces drop below the cutoff, just terminate the
run. For frequencies only the residual forces matter.
And on this note, if you need to get frequencies for
constrained optimized structure, I have a small program
for level-shifting the frozen atoms away before
determining the frequencies. That way you get 'real'
frequencies within the constraints that you have put on the
system, i.e. all contributions from frozen atoms are
completely decoupled from the remaining degrees of
freedom. If you are interested drop me a note.

        Frank

--------------------------------------------------
| Frank Jensen, Department of Chemistry          |
| SDU Odense University, DK-5230 Odense, Denmark |
| FAX +45 66 15 87 80 , Voice +45 65 50 25 07    |
| http://bogense.chem.ou.dk/~frj                 |
--------------------------------------------------

On Thu, 30 Nov 2000, Gerald Lushington wrote:

> Hi again!
> 
> No intention of closing down this discussion so quickly if others have points
> to make, but I thought I'd take a moment and share the three very insightful 
> responses I've gotten so far.  Each response addressed different aspects of the
> problem and in each case proved very helpful.  Glad to see you're all so 
> coordinated  ;)
> 
> Jeremy Greenwood confirmed that convergence is slow, but also that it's 
> generally possible to (eventually) achieve.  You may wish to refer below to his
> response for a few constraining tips he suggests.
> 
> Nicolas Ferre gave a few suggestions as to how to circumvent the strange 
> "Delta X" and "New X" values we were getting.  Turns out that using the "NOSYM"
> keyword straightens this glitch out nicely!
> 
> Finally, as you have likely seen, Stefan Fau provided a very thoughtful and 
> detailed discussion on the validity of frequency analysis on our constrained
> system.  Have to sit down and think carefully about this, but his suggestions
> make sense to me on the first pass.
> 
> Anyway, thanks all!  Original question and the three responses appended below.
> 
>                                   - Gerry
> 
> ----------------------------------------------------------------------
> 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
> 
> 
> 
> == Original Question ==========================================================
> 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.
> 
> 
> == Response #1:  Jeremy Greenwood <[email-masked]> ====================
> 
> I am doing the same kind of thing at the moment.
> 
> I *have* been able to get the displacement all the way down eventually, 
> though it gets harder the larger the system and the greater the constraints.
> 
> Are you constraining cartesian co-ordinates? I suggest trying to constrain
> a minimal number of internal co-ordinates for key heavy atoms.
> 
> Hope this helps a little,
> 
> Jeremy
> ----------------------------------------------------------------------
> Jeremy Greenwood                                 [email-masked]
> Department of Medicinal Chemistry                      bh +45 35306117
> Royal Danish School of Pharmacy                        fx +45 35306040
> Universitetsparken 2, DK-2100 Copenhagen, Denmark      ah +45 32598030
> ----------------------------------------------------------------------
> 
> 
> == Response #2:  Nicolas Ferre <[email-masked]> ================
> Hi Gerald,
> 
> I'm very surprised by the output you give. Normally the total Delta X is
> the sum of the 2 previous one, Linear and Quad. And when you freeze some
> variables, all of theses 3 quantities have to be zero. I've performed
> some calculations on crambin, allowing only 3 residues to move, so
> freezing other 43 residues (it's a QM/MM calculation) and the output
> looked fine.
> In your case, it seems to me that only R1 to R5 are really frozen, but
> that atoms 1, 11, 38, 48 and 50 are translated by the vector (2.501,
> -8.90656,8.10547) in bohrs. I don't understand why.
> Are you sure that these atoms don't move ? Did you check the values
> during optimization ?
> How did you freeze the coordinates of these atoms ? With a "-1" in the
> How did you freeze the coordinates of these atoms ? With a "-1" in the
> Z-Matrix or with the MODREDUNDANT keyword ? Did you use the NOSYM keyword
> (optimization in the input orientation or in the standard orientation) ?
> 
>                                         Nicolas
> 
> 
> == Response #3:  Stefan Fau <[email-masked]> =================================
> 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
> 
> ===============================================================================
> 
> 
> 
> -= This is automatically added to each message by mailing script =-
> [email-masked] -- To Everybody    |   [email-masked] -- To Admins
> [email-masked] -- HELP CHEMISTRY or HELP SEARCH
> [email-masked] -- archive search    |    Gopher: gopher.ccl.net 70
> Ftp: ftp.ccl.net  |  WWW: http://www.ccl.net/chemistry/   | Jan: [email-masked]
> 
> 
> 
> 
> 

----------------------------------------------------------------------
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