RE: AMBER: more than 100,000 atoms ?

From: Kristina Furse <>
Date: Fri, 8 Aug 2003 13:14:10 -0500


I think your choice of method here depends pretty heavily on what kind of PMF
you want to do. In my case, I want to bias the sampling for a particular
torsional angle of my substrate inside my enzyme. I, too, checked out Gibbs,
first but really wanted (more like needed with a system this big) the PME
improvements that came with amber7. Since I'm driving a torsional angle, I
could just use NMR-style constraints in sander to bias my sampling, so I guess
if your PMF involves biasing a distance, angle or torsion this will work. But
if it's something else, you may need to find a way to use Gibbs.

There are one or two nice emails in the reflector archive (esp. one from Case
on Fri Aug 16,2002 Re:PMF with sander module, and from Pratul Agarwal on Wed
Aug 14, 2002, same subject line--I have printouts, or I'd send the actual
link), that discuss PMF in sander. They have references for the introduction
of the WHAM method of removing the bias to generate the PMF curve which were
essential for background on the method (Kumar, et al. J. Computat. Chem, 1992,
13:1011-1021 and Kumar et. al., J Computat Chem, 1995, 16:1339-1350), but
didn't lead me to any executable I could use (I'm not at the point where I can
write my own...). A google search lead me to Alan Grossfield's website
( where I found a nice WHAM executable. I tested
it on sander output (Alan uses another platform) where I drove the ethane
rotor, and it gave very reasonable results considering I ran a pretty
rudimentary test of ethane with gB solvent so I couldn't really control the
"temperature." Most important lesson learned the hard way: amber expects
torsional constraints in kcal/radian**2, and Grossfield's program expects
constraints in kcal/degree**2, so a little conversion is needed.

Anyway, if you think that you can and do want to use sander7 or pmemd to run
you PMF, let me know and I can tell you more about what I did.

Good luck!


>===== Original Message From David Smith <>
>Hi all,
>Thanks to Tom and Robert for their quick responses.The reason I have
>specifically CC'd this to Kristina will become clear shortly.
>Actually, I was using tleap and (apparently) I had no problems with this
>stage. That is, both the prmtop and crd files have an I6 natom. I didn't
>sift through the code to find out why.
>The real problem was sander6. I had a look at sander7 (should've done it
>before, I know) and now understand (a bit) how it is handled there.
>My main problem is that I want to do some PMF calculations on my large
>system. I had previously been using GIBBS for these.
>I note that gibbs7 still reads I5. I tried to change gibbs6 to read (and
>write) I6. That is, I changed the 9018 format line in GETCOR
>and MDWRIT in giba.f to FORMAT(I6,5E15.7). This allowed me to start
>running (no complaints about number of atoms not matching between prmtop
>and crd files) but then I just get zero's for everything (all coords and
>energies etc). I admit my programming is not that great so I am a bit
>lost as to how to fix that.
>If I would like to do my PMF's on the big system, I basically see three
>(a) get gibbs6 (or 7) running to read I6.
>(b) Use PMEMD in a way similar to what Kristina was doing.
>(c) Use Sander 7 in a way similar to what Kristina was doing
>The first problem with options b and c is that I don't have any
>experience with doing PMFs this way. In this respect, I was wondering if
>Kristina could give me a brief description (maybe an example) as to how
>it is done and what the WHAM method is all about. Maybe a reference or
>something that could lead me in the right direction.
>I guess the second problem is whether or not the restraints are working
>properly in PMEMD. Therefore, I would really appreciate it if any
>developments by Robert and Kristina in this direction could be posted to
>the list.
>Finally, if anybody has any suggestions as to how to get gibbs running
>with I6, I would also immensely appreciate that.
>I think I have now spent my adjective database and will politely sign
>Thanks to everyone for their patience.
>On Thu, 2003-08-07 at 22:33, Robert Duke wrote:
>> David -
>> Sander 7 can do this, no problems. So can PMEMD, if you are mostly still
>> an amber 6 context. The prmtop is not a problem but your first inpcrd file
>> is. Really all that is needed is the change of natom from i5 to i6; if you
>> are coming from tleap or xleap, I don't know how it normally handles this.
>> I whacked the xleap source to get it to do the right thing, myself, but
>> is not a standard part of the distribution. The change I made is shown by
>> the following
>> diff (amber 6 source):
>> diff -tr ./leap/unitio.c /home/rduke/amber6/src.ori/leap/src/leap/unitio.c
>> 5466c5466,5478
>> < FortranFormat( 1, "%5d" );
>> ---
>> > /*
>> > * Write a amber 6 inpcrd atoms line for less than 100,000 atoms, and a
>> > * amber 7 inpcrd atoms line for 100,000 or more atoms. This allows us
>> > * use PMEMD with more than 99,999 atoms.
>> > */
>> >
>> >
>> >
>> > if (iVarArrayElementCount(uUnit->vaAtoms) < 100000)
>> > FortranFormat( 1, "%5d" );
>> > else
>> > FortranFormat( 1, "%6d" );
>> >
>> The guys I work with routinely deal with systems of over 100,000 atoms,
>> which is one reason I modified sander 6 into pmemd (faster, handles bigger
>> problems (up to about 329844 atoms)in less memory). See the amber web
>>, for the link to get a copy of PMEMD.
>> Regards - Bob Duke
>> ----- Original Message -----
>> From: "David Smith" <>
>> To: <>
>> Sent: Thursday, August 07, 2003 3:44 PM
>> Subject: AMBER: more than 100,000 atoms ?
>> > Hi all,
>> >
>> > Apologies if this question has come up before but I don't remember it
>> > and a quick search of the archives didn't net me anything.
>> >
>> > I wanted to try a system with more than 100,000 atoms. I know it's a bit
>> > big but my boss is an ambitious guy.....
>> >
>> > I gave it a go and initially got very confused by getting caught by
>> > "Peek Ewald" in a simple sander minimization.
>> >
>> > After a bit of surprise and a few write statements I realized that NATOM
>> > was being read as "I5" from the CRD file. Changing it to "I6" got me
>> > past "Peek Ewald" but I then got stuck a bit later. The web page tells
>> > me that, in general, NATOM is "I6" in PARMTOP but "I5" in CRD.
>> >
>> > I am still "trapped" in Amber 6 but my guess is that these formats
>> > persist into Amber 7 ?
>> >
>> > Can anyone recommend a way to do such a system with AMBER (any version)?
>> >
>> > i.e. Would I have to change all the formats associated with "read (9,F)"
>> > statements. Is there a semi fast way to do it ? Even if I did manage
>> > that, would the developers imagine that something else would break and
>> > that I should simply abandon the whole idea.
>> >
>> > Thanks for your patience ....
>> >
>> > David.
>> >
>> >
>> > --
>> > ---------------------------------------
>> > Dr. David Smith
>> > Department of Chemistry
>> > Ludwig Maximilians University
>> > Butenandt-Str. 5-13, D-81377 Munich
>> > Germany
>> > Tel.: +49 (0)89 2180 77740
>> > Fax.: +49 (0)89 2180 77738
>> > e-mail:
>> > ---------------------------------------
>> >
>> >
>> > -----------------------------------------------------------------------
>> > The AMBER Mail Reflector
>> > To post, send mail to
>> > To unsubscribe, send "unsubscribe amber" to
>> >
>> >
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to
>> To unsubscribe, send "unsubscribe amber" to
>Dr. David Smith
>Department of Chemistry
>Ludwig Maximilians University
>Butenandt-Str. 5-13, D-81377 Munich
>Tel.: +49 (0)89 2180 77740
>Fax.: +49 (0)89 2180 77738

Kristina E. Furse
Department of Chemistry
Center for Structural Biology
Vanderbilt University

The AMBER Mail Reflector
To post, send mail to
To unsubscribe, send "unsubscribe amber" to
Received on Fri Aug 08 2003 - 19:53:01 PDT
Custom Search