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.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sat Aug 16 2003 - 21:53:02 PDT
