Implementation of Gay-Berne ellipsoids

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I'm trying to understand how Gay-Berne ellipsoids are implemented in LAMMPS. I've read the documentation at http://lammps.sandia.gov/doc/PDF/pair_gayberne_extra.pdf, but when I do calculations, the results don't match what I'm expecting based on that document. So I started going through the source code (pair_gayberne.cpp) and comparing to the equations in the documentation, and it really looks to me like they don't match. Perhaps someone can explain to me what's going on?

The first difference is in the definition of the B matrix. According to the documentation, B = A^T E^2 A, where E = diag(epsilon_a, epsilon_b, epsilon_c). But here's how it's computed in the code (in PairGayBerne::compute()):

      MathExtra::diag_times3(well[itype],a1,temp);
      MathExtra::transpose_times3(a1,temp,b1);

"well" is computed in PairGayBerne::coeff() as:

        well[i][0] = pow(eia_one,-1.0/mu);
        well[i][1] = pow(eib_one,-1.0/mu);
        well[i][2] = pow(eic_one,-1.0/mu);

So instead of E^2 it's using E^(-1/mu).

I also found differences in how it computes dU/dq and dChi/dq. According to the documentation each of these involves a multiplication by A, but the code (in PairGayBerne::gayberne_analytic()) seems to be leaving that out.

Thanks in advance for help in understanding this!

Peter

2

Mike Brown wrote that pair style, and can likely answer.

Steve

3

Hi Peter,

Gay-Berne is a tough cookie. I remember that it was painful to get the code implemented correctly into LAMMPS (I actually did a dual implementation in Matlab and LAMMPS to help debug) and have helped several on this mailing list with reproducing. It has been tested by the group who developed the biaxial formulation and compared to their code and I am not aware of any bugs in the implementation.

The relationship between the variables in the code and the documentation should be treated very loosely because it is not optimal to do a literal implementation of the math sequences as described. For your first question, for example, see the definition of chi_12.

Best of luck!

  • Mike
4

Hi Mike,

Chi_12 seems to be implemented exactly as described in the document:

  tempv[0] = iota[0]/r;
  tempv[1] = iota[1]/r;
  tempv[2] = iota[2]/r;
  double chi = MathExtra::dot3(r12hat,tempv);
  chi = pow(chi*2.0,mu);

Consider this question: if I increase eps_a, eps_b, and eps_c, should that increase or decrease the strength of the interaction? Unless I'm totally missing something, the equations in the document clearly show it should decrease the strength. Increasing the epsilons means increasing E, which increases B, which decreases Chi. But in the code it's implemented the opposite way: increasing the epsilons *decreases* E, and therefore *increases* the strength of the interaction. And I've verified that it really does have that effect.

So it's not just a question of the math being rearranged. The implementation and the documentation disagree about the effect of those parameters.

Peter

5

K, I see the problem. The LAMMPS behavior is the desired behavior. The gayberne_extra documentation is inconsistent. The epsilon_i_. parameters specified in LAMMPS are the well-depths for the side-to-side, end-to-end, etc. interactions (not inversely proportional as listed in the extra documentation) such that

E_i=diag{e_i_a^(1/mu), e_i_b^(1/mu), e_i_c^(1/mu)}

I will fix the document. also note that the Everaers reference on the LAMMPS gay-berne documentation page is probably best for understanding both the Gay-Berne and RE-squared potentials (and I tried to follow this notation)

Thanks for reporting. – Mike

6

Thanks! That makes sense. I guessed it was probably something like that.

While you're fixing the document, do check the equations for dU/dq and dChi/dq as well. I think each of them includes an extra multiplication by A that isn't in the code.

Peter