Subject: AMBER: frozen x unfrozen and the gibbs tragedy
From: David Smith
Date: 23 Oct 2003 17:56:00 +0200

Hi Herbert,

This may be of little help but the gibbs program has the possibility to
add holonomic (SHAKE-type) constraints to arbitrary internal
coordinates. In this sense it is already in the code.

To use this facility for a dynamics run you would need to:

(a) define a "null" perturbation topology file (saveamberparmpert in
leap with perturbations set to false in the unit definition, p192 of
amber 7 manual).

(b) set NSTLIM to the number of steps you want to do, but make NSTMEQ
bigger than this number (p160).

(c) set intr=1 in the gibin file.

(d) use the "formatted" input to define your constraints (p184) e.g.:

00101 00102 00000 00000 00000 00000 00002 000.00000 001.00000
0000.00000 003.50000 000.00000 003.50000 0000 0000

would constrain the bond between atoms 101 and 102 to be 3.5 Angstroms
(if you set almda=1.0 and isldyn=-3 you only really need the second 3.5
but why not have both for emphasis). You can also consider 999 in place
of 3.5 which will then use the value found in the PINCRD

(e) run gibbs to your hearts content.

Warnings.

(1) You should read and understand p186 if you want to constrain
dihedrals.

(2) If you have lots of constraints, you may want to write a little
program to prepare the input for you. This is pretty trivial as all the
formats you will need are given on p184.

(3) It seems that almost nobody uses gibbs anymore. Certainly nobody is
looking after the development. This means that, in particular if you are
running a periodic system with ewald on, you will have some trouble.
That is, ewald is turned off by default in amber 7. If you want to use
it you can recompile it with FLAGS = -DEWALD in the gibbs Makefile. If
this works for you then good, but you can forget parallelization because
the combination of ewald and mpi is disabled. Alternatively you can use
gibbs from amber 6 in which ewald is enabled but mpi is globally
disabled.

(4) In the case that you are really persistent and do go with one of
these options, there is just one more additional warning. As the ewald
implementation is an old one (e.g. you will need to use some funny
formats for the ewald input which are not even documented in amber 7
manual but they are in amber 6 for both gibbs and sander_classic), you
will want to "re-equilibrate" your system in gibbs before you collect
too much data. This is especially true if you are running NPT but is
probably also necessary for NVT.

So, if you are feeling very brave go ahead with the gibbs approach.

In my humble opinion, the above situation is another example of what I
have come to think of as the great gibbs tragedy. That is, there are
many such great features in gibbs which are virtually unusable these
days because the dynamics code didn't keep up with the advances in
sander and pmemd. I fully understand that it is a huge effort to bring
it into line and I certainly can't put my own hand up to do it. Hence,
it seems as if gibbs is destined to fade away, underused and
underappreciated (please excuse the melodrama, I'm just a bit bored
today).

My two cents.

David.



On Thu, 2003-10-23 at 16:25, Herbert Georg wrote:
> >
> >
> >>I would like to perform a simulation with a system in which one molecule
> >>has to be flexible and the other has to be frozen (at least the bonds
> >>have to be frozen). Is this possible??
> >>I was reading about the keyword NTR and the group section but it's not
> >>clear to me if the restraints are about the absolute positions of the
> >>specified atoms or if its about their bonds.
> >>
> >>
> >
> >Ntr refers to the abosolute positions. I'm not sure what you mean when you
> >say one molecule has to be "frozen". You could use ntr=1 to keep one molecule
> >in a (nearly) fixed position and allow the other one to move.
> >
> >..good luck...dac
> >
> >
> By frozen I mean that the molecule must have not internal degrees of
> freedom (at least for bond lengths). While the other molecule must have
> internal degrees of freedom.
>
> >
> >
> >>Can I apply SHAKE to the molecule I want to keep frozen only?
> >>
> >>
> >
> >You can apply SHAKE to the part of the molecule that is also restrained by ntr
> >restraints.
> >
> >..good luck...dac
> >
> >
> >
>