AMBER Archive (2009)

Subject: [AMBER] igb=0 vs. igb=6 give different answers

From: Hugh Heldenbrand (helde010_at_umn.edu)
Date: Thu Apr 09 2009 - 15:51:02 CDT

This problem is related to a previous thread
(http://structbio.vanderbilt.edu/archives/amber-archive/2009/1641.php),
but I thought that I would start a new one since it could be a different
issue.

I have been doing gas phase minimizations of a small residue of my own
creation and a TIP4PEW water molecule. To do a gas phase minimization I
usually set igb=0, but I often have failures due to my system exceeding
the limits of the virtual box (see previous thread). Dr. Case suggested
that using igb=6 should avoid this problem. In the AMBER manual it says
that these two settings should be logically equivalent and there is a
previous listserv post about this as well
(http://structbio.vanderbilt.edu/archives/amber-archive/2008/2647.php).

However, I see a 10 kcal difference in energy when I minimize my system
with igb=0 vs 6.

My input file is pretty straightforward:

Energy of the A-H2O interaction
&cntrl
  imin = 1, !Perform minimization and no MD
  ntpr = 1, !Print every step
  maxcyc = 50000, !Number of minization steps (no, it never takes
this long to reach the convergence criterion)
  ntb = 0, !No periodic boundary
  igb = 0, !No implicit solvent
  cut = 999 !No cutoff
/

I looked carefully at the other default flags that we being used in my
two output files. I noticed that when igb=6 there were several born
solvation flags popping up that weren't there when igb=0, such as:

     saltcon = 0.00000, offset = 0.09000, gbalpha= 1.00000
     gbbeta = 0.00000, gbgamma = 0.00000, surften = 0.00500
     rdt = 0.00000, rgbmax = 25.00000 extdiel = 78.50000
     alpb = 0

But those don't seem to be actually involved in the calculation (I tried
setting surften = 0.0000 and extdiel = 1.0000 and my answer was still
the same and still 10 kcal different from when igb=0).

The other difference that I noticed was that when igb=0, these options
that seem to be related to the extra point of my TIP4PEW would be set
automatically: frameon=1, chngmask=1, followed by a list of different
terms being trimmed:

| EXTRA_PTS, trim_bonds: num bonds BEFORE trim = 10 0
| EXTRA_PTS, trim_bonds: num bonds AFTER trim = 10 0
| EXTRA_PTS, trim_bonds: num bonds BEFORE trim = 14 0
| EXTRA_PTS, trim_bonds: num bonds AFTER trim = 13 0
| EXTRA_PTS, trim_theta: num angle BEFORE trim = 13 0
| EXTRA_PTS, trim_theta: num angle AFTER trim = 13 0
| EXTRA_PTS, trim_theta: num angle BEFORE trim = 19 0
| EXTRA_PTS, trim_theta: num angle AFTER trim = 19 0
| EXTRA_PTS, trim_phi: num diheds BEFORE trim = 20 0
| EXTRA_PTS, trim_phi: num diheds AFTER trim = 20 0
| EXTRA_PTS, trim_phi: num diheds BEFORE trim = 28 0
| EXTRA_PTS, trim_phi: num diheds AFTER trim = 28 0

However, in the igb=6 case frameon =0, chngmask=0 are the defaults. I
tried to create an ewald namelist to change those options to 1, but it
seems that they were ignored, the calculation still uses frameon =0 and
chngmask=0 and returns and energy that is 10 kcal lower that igb=0. Is
igb=6 treating my extra point as a real atom bonded to the oxygen
water? Or is something else going on?

Thanks,
-Hugh Heldenbrand
Graduate Student
Chemistry Dept.
University of Minnesota

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber