Trouble with lj units

LAMMPS LAMMPS Beginners
1

Dear LAMMPS Users,
The thing is I’m having trouble understanding how lj unit system and the interconversion of real to lj and lj to real units works. I thought I knew. I went through its documentation several times but I’m still having hard time figuring out very basic stuff on my own. So kindly validate whether I’m going in right direction or not.

Just like in SI system we have 7 fundamental quantities with mass, distance, time etc. and everything else is derived from these. Similarly, In lj unit system, we have m, \sigma and \epsilon and rest other quantities are derived from these three. mass is written in m units, distance in \sigma units and energy in \epsilon units. By default LAMMPS sets the value of m, \sigma and \epsilon (& k_B as well) fundamental quantities to 1. So far so good?

Now, Suppose I have a graphene sheet with 100 atoms. The molar mass of carbon is 12.011 g \space mol^{-1}. So, I’ll simply calculate per atom mass of carbon which is 1.99 e-26 Kg but since we have to make quantities dimensionless for lj unit system so we have to divide it by m (which would also be in kg)? This is where my understanding gets a jolt. Just to be clear, should I take M to 100*1.99 e-26 Kg and set m to 1.99 e-26 Kg which makes my M^* to be 100. Am I correct in making this assertion? (Since the formula is M^* = \frac{M}{m}). Also, what if I have multiple atoms of different masses?

Moving on, at equilibrium distance r = 2^{1/6} \sigma so based on interatomic distances in the graphene I can set my \sigma = \frac{r}{2^{1/6}} accordingly. Through this I can define my 1*\sigma and based on this rest all the distances will defined in the units of \sigma. Right?

The corresponding Energy value at this, r = 2^{1/6} \sigma will be my epsilon?

Lastly, since the time unit depends on m, \sigma and \epsilon. Suppose I want to run my simulation for 200ps. So I’ll convert it into respective lj units by \tau^* = \frac{200*10^{-12}}{\sqrt{\frac{\epsilon}{m\sigma^2}}}.

P.S. I understand @akohlmey made it very clear here to not do MD on your own but since I’ve already started this venture. I can’t not learn just because I don’t have anyone to learn from.

2

Hello,

For a graphene system, your life would be so much easier if you were using the unit style real or metal. LJ units are rarely a good choice in those situations… Is there a reason for you to use LJ here?

I wrote a tutorial on graphene, may be you could be interested.

Best,
Simon

1 Like
3

What could be done instead is to use the mass of a carbon atom as the reference, so m = 1 (in LJ units) would corresponds to the carbon atom. Then, if you have another type of atom in your simulation, let us say a hydrogen atom, its mass in LJ units would be 1.008/12.011=0.0839.

1 Like
4

Thanks a lot @simongravelle. I’ll surely check it out. I’m not simulating graphene. I used it as an example to explain where I lack conceptual understanding.

5

Can you please shed some light on \epsilon & unit of time in lj units as well?

1 Like
6

Its the exact same principle. If you set carbon as the reference, then \epsilon_\text{CC} = 1 and \sigma_\text{CC} = 1.

1 Like
7

Hi @Akshay_Chauhan,

I step upon your thread and though that I could give one insight or two to your situation.

  • Reduced units are a very common practice in physics (more than chemistry), and the scope of your problem ranges outside MD. So I think your question is still relevant for anyone struggling to wrap their head around it.
  • Maybe taking the problem in reverse might help: from a set of lj units, how can you go back to some physical quantities? As @simongravelle pointed out, the idea is to take a reference mass/time/energy and express all masses/times/energies with regard to your reference. That is, from the real units set, notice how we can directly divide g\cdot{}mol^{-1} quantities. The important part is only the ratio of one with regard to your reference (in a given set of units). That’s all. The same goes for other reference quantities. The choice is usually made from quantities easily derived from the model you use (bond length, force constant, etc…). It depends mostly on the physical properties you are interested in ans the model you are using.
  • That means if you take some m, \varepsilon and \sigma (from a given model giving you the values) as your references, your reduced time unit in LAMMPS will be \tau = \sqrt{\frac{m\sigma{}^{2}}{\varepsilon{}}}. It is up to you to compute this value in a given unit set when you want to convert your results back. (Note that the default time step is 5\cdot{}10^{-3}\tau so unless you change it, \tau is reached after 200 integration steps.)
  • I read @akohlmey 's advice you referred to a bit differently than “do not do MD on your own” in the sense that “do not learn if you don’t have anyone to learn from”. It depends where you stand and what you expect from learning MD. No one can stop you to learn on your own and ask questions if you’re interested in the underlying physics/programming. You can easily play around with LAMMPS, make a double pendulum or classical planetary systems. But if you expect professional research outputs from MD simulations, you will need personal guidance, and I agree with that. In all cases you have to be aware that guidance from forums and mailing lists (even though they contain a lot of info) is limited, and that the best way to learn will always be personal interaction. At some point, if you want to learn “MD craftsmanship”, you’ll have to actively look for some “craftsman” accepting to teach you.
4 Likes
8

Dear Lammps Users,
I’m having slight trouble with understanding dielectric value for my lammps script. Based on the reference \sigma and reference \epsilon (energy), the value of 1e (e being electronic charge) comes to be 15 in LJ units. The paper I’m currently trying to reproduce have set their Bjerrum Length to be 3 \sigma. Based on the definition of Bjerrum Length, \lambda_B = e^2/(4 \pi \epislon_0 \epsilon k_B T). (\epsilon is the dielectric in this case)

The formula in reduced units becomes, \lambda_B = q^2/(\epsilon). Using which the dielectric value comes to be 75. According to the documentation, " a value of 4.0 reduces the Coulombic interactions to 25% of their default strength" a value of 75 seems very unphysical because the coulombic interaction reduces to 1.3% of their default strength.

Is my line of reasoning correct?

Yours Sincerely
Akshay

P.S. I didn’t want to start a new thread for this, hence using this old thread since it too pertains to lj units.

9

Hi @Akshay_Chauhan,

A few comments I can make:

  1. You can make your formula more readable by enclosing your LaTeX writing in $ symbols. There is a LaTeX interpreter on the forum
  2. The Bjerrum length formula is: \lambda_B = e^2/(4 \pi \epsilon_0 \epsilon_r k_B T)., notice that it does not depend on the units you use. You could compute it using some units and divide it using your unit scale to get the LJ value, that’s the easiest way. However I think your reasoning is incorrect, k_B is set to 1 in LJ units, but T and \epsilon_0 also have values to scale, and \epsilon_r (the dielectric constant) depends on your medium. But that’s not the parameter I would try to modify first as electrostatic interactions are generally simulated in vacuum, with \epsilon_r=1. Other values are for implicit solvent and other specific cases. With this in mind, I think your formula in reduced units is incorrect (and non-homogeneous to a length).
  3. What do you mean when you say “The paper I’m currently trying to reproduce have set their Bjerrum Length to be 3 \sigma”? What are you trying to simulate? This is not something that we can guess from the informations you provide.
  4. You mention some part of the dielectric command manual page. Please state when you do so. Your statement is a bit confusing.
  5. This is unrelated to your previous issues with LJ units. You should start a more detailed thread next time stating exactly your specific problem.
10

My apologies for the ambiguity in stating my problem and all other confusions. I have created a new post for the aforementioned problem.

Thanks and Regards
Akshay