AMBER Archive (2003)

Subject: AMBER: minimum image convention mystery

From: Martin Mucha (mmuc7376_at_mbox.troja.mff.cuni.cz)
Date: Sat Aug 16 2003 - 15:42:37 CDT

Dear Amber users,

I have been trying to simulate binary mixture of Lennard-Jones
particles in a periodic box, in the NVE ensamble. However, the energy
did not conserve, but went wild instead. To be more specific,
suddenly, the kinetic energy increased, as if a result of large
forces. I am sure that my system was equilibrated, actually
I took an equilibrated coordinates from a gromacs run. Therefore
I put a few lines into the code of 'subroutine short_ene', which
is basically the inner loop of LJ forces calculations, that would
report if a LJ force between two particular particles was greater
than, say, 10 (of whatever the internal units of force amber uses),
alongside with the indices of the two bad-behaving atoms in question.
Then I looked at the first occurence of an unphysically large force,
and was further investigating what is the cause. Of course,
the strightforward reason is that the two atoms in question
got very closely together(being only 1 angstroem apart). However,
the question is how did they get that close. I found out,
that while they were approaching each other, the LJ force was zero,
and was suddenly "switched on" only when they were so close, that
the repelling force was enormous. It seems to me, that it has something
to do with the nonbonded neighboursearching or minimal image convention
calculation. Let me give you some more evidence:

In the inner loop, more precisely in the file "ew_directp.h",
I have inserted these three lines of code, right after the
square of the distance "delr2" is calculated

C------------------------------------------------------------
      if ((i.eq.AT1.and.j.eq.AT2).or.(i.eq.AT2.and.j.eq.AT1)) then
         write (44,*) 'inner ',sqrt(delr2)
      endif
C------------------------------------------------------------

Note that AT1 and AT2 should be replaced by the indices of the
two bad-behaving atoms I am investigating, let say:
#define AT1 508
#define AT2 25

Using this code I can find out, what does the inner loop think
is the distance between the closest images of the
two atoms.

I compared it to the distance between the two atoms, as calculated
directly from the coordinates, I put the folowing lines of code
into subroutine runmd in "runmd.f", into loop over MD steps,
right before the force caluclation (CALL FORCE(...)):

C--------------------------------------------------
      write (44,*) 'STEP ',NSTEP
      write (44,*) 'outer ',sqrt((X(AT1*3-2)-X(AT2*3-2))**2+
     +(X(AT1*3-1)-X(AT2*3-1))**2+(X(AT1*3)-X(AT2*3))**2)
C--------------------------------------------------

Of course, if the two atoms, not just their images, are close
to each other, then the "outer" and the "inner" distances should
be equal. However, they are not, look at the output from the code
above: 

---
 STEP  0
 outer   4.00686354
 inner   4.00686354
 STEP  1
 outer   3.9797466
 inner   3.9797466
 STEP  2
 outer   3.95117027
 inner   3.95117027
 STEP  3
 outer   3.92125795
 inner   3.92125795
 STEP  4
 outer   3.89019832
 inner   3.89019832
 STEP  5
 outer   3.85802023
 inner   44.0021111
 STEP  6
 outer   3.82411395
 inner   43.9659115
 STEP  7
 outer   3.78857943
 inner   60.2610321
 STEP  8
 outer   3.75140294
 inner   56.2451189
 STEP  9
 outer   3.71270959
 inner   56.252646
---
Until the 4th step they are equal, but on the 5th step
the inner loop thinks that the two atoms are very far apart (44.0021111),
although actually they are very close (outer   3.85802023).
The wrong distance calculated within the inner loop is caused
by taking the wrong images of the two atoms.
  We can easily guess what happens next, the two atoms just continue
to fly approaching each other, because they do not feel the reppeling
LJ force, since the inner loop thinks that they are far apart.
After a while the inner loop changes it's mind and finally 
realizes that they,
are close to each other, but now they are too close already,
and the force is enormous, they get a kick, eventually warming up,
the whole system, the energy goes wild.
  So my question is: Is there something wrong with the construction
of neighbour list from which the minimal images are determined? A bug in
Amber?

For completeness:
 I am using Amber 7, the same happens in Amber 6.
 time step is pretty reasonable for a LJ system, dt=0.01
 LJ parameters are those for Argon and Krypton.
I calculate the neighbour list every step:
    nsnb=1
no pme or other tricks:
   ew_type=0
   vdwmeth=0
   eedmeth=1
   use_pme=0
   nbflag=0
wrapping the coordinates:
   iwrap=1
(setting iwrap=0 will not change the problem)

   Martin Mucha

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu