[lammps-users] meam: derivative of rhobar1(and 2) for mixed pair potential

LAMMPS LAMMPS Mailing List Mirror
1

Steve or Greg,

I’ve looked through the MEAM code and am unable to find where it computes d (rhobar1)/dr. This is the background density in a mixed pair setting. Currently in LAMMPS it is a special case for the C11b lattice type. Is there some sort of interpolating going on or am I just missing it?

Also, are there any plans to include the L12 crystal structure in the library?

Thanks,
Michael Sellers

2

Dear all:

    In my simulation, I found an inconsistency between the early version of
Lammps (July 16, 2008) and the newer version (Jan 19, 2009) regarding
temperature computation.
    My system is made of "bcb" and "ch2" molecules, where "bcb" molecules are
rigid bodies (basically phenyl) and "ch2" molecules are particles(CH2 united
atom).
    I use the command below to do the thermostat:
    fix 1 bcb temp/rescale 10 15 15 0.05 1.0
    fix 2 bcb rigid molecule
    fix 3 ch2 nvt 15 15 1 drag 1.0

    The early version of Lammps gave the thermo temperature of the system to
be about 15, which seems to be reasonable because I set the equilibrium
temperature to be 15 for both groups. However, the latest version computed the
thermo temperature of the system to be about 22. And, when I used "compute
temp" to compute temperature for "bcb" and "ch2" respectively, I found
temperature for "bcb" was about 15 while temperature for "ch2" was about 12.
    I was confused by the result. Could anybody tell me what causes this
problem?

    Thanks!

jiwu

3

Greg will have to answer these Qs.

Steve

4

Computing the temp and thermostatting rigid bodies
is tricky. So I'd focus on the current version and whether
it is right. (the past version may have had its own problems,
and things have changed a bit with temp/thermostatting
since then).

In general, I'd use fix langevin for rigid bodies, not temp/rescale.
But first I would run the whole system NVE and see if the
2 components equilibrate to the same temp. And I'd viz
the system to make sure it is doing what you expect.
If this fails, then please make a small system (e.g. 50 molecules
of each in a box) and post it and I'll take a look.

Steve

5

Thanks for your advices.
  I ran a simulation without thermostatting, namely, NVE for both groups. They
failed to equilibrate at the same temperature.
  After checking the code in fix_rigid.cpp, I ' m afraid it has to do with
this change in the latest lammps:

1105 int n = 0;
1106 if (domain->dimension == 3) {
1107 for (int ibody = 0; ibody < nbody; ibody++) {
1108 n += 3*nall[ibody] - 6;
1109 if (nall[ibody] == 2) n++;
1110 }
1111 } else if (domain->dimension == 2) {
1112 for (int ibody = 0; ibody < nbody; ibody++)
1113 n += 2*nall[ibody] - 3;
1114 }

    At the line 1108, the old code was:
    if (nall[ibody] > 2 ) n+= 3*nall[ibody] - 6;

    I tried to print out "nall[ibody]". Sometimes the value is zero ( probably
for the group that is not rigid body ?). After I added the condition " if
(nall[ibody] > 2) " back, the system can equilibrate again.
   
    Also, there is a typo in the line 230:
230 for (int m = mlo; m <= mhi; i++) {

      where i++ should be m++. In the case when keyword "force" is used with
fix_rigid, lammps enters infinite loop due to this bug.

jiwu

6

Yes - I think you nailed the bug. The code was
not recognizing that the temperature group might not
include any atoms in some rigid bodies. The change
to handle 3d vs 2d and dimers correctly introduced this bug.

I posted a patch - 22Jan09 - see if that now works
for your system. Also see the enhanced discussion of
DOF on the fix rigid doc page.

Thanks,
Steve