[lammps-users] Calculation of heat flux

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,

I want to calculate the heat flux in a solid system with periodic boundary conditions. The microscopic formula I’m implementing in a compute is as follows:

J = sum_i v_i epsilon_i + 1/2 sum_pairs F_ij r_ij . (v_i+v_j)
where:
v is velocity
epsilon is local energy (KE and an appropriate share of PE)
F_ij is the force that atom j exerts on atom i
and sum_pairs means that i and j are iterated such that each pair is visited once only.

The neighbor lists in LAMMPS are (at least for pairwise interactions) constructed such that a double sum over each real atom i and all their neighbors j visits each pair once, so that aspect is fine. However, I have a problem if j happens to be a ghost atom, i.e. a periodic copy of a real atom, since LAMMPS only keeps track of the positions of ghost atoms but not their velocities.

How can I find out the velocity of a ghost atom? Two posssibilities occur to me, but I can’t make progress with either:

  1. I assume that when positions are being updated after an integration step the positions of real atoms are simply copied to their corresponding ghosts, i.e. if there are 30 real atoms and atom 43 is a (periodic) copy of atom 7, I would expect there to be a line of code like x[43][0] = x[7][0]. I’ve looked but can’t find it. Where is it and can I simply copy this book-keeping for the velocities, writing v[43][0] = v[7][0]?

  2. Somewhere there must be a list recording that 43 is simply a copy of 7, i.e. something of the form parent[43] = 7 (which would also contain the information parent[7] = 7). When in my formula I need v[43], how can I access this list to perform the necessary mapping? Instead of v[j][0] I could then write v[parent[j]][0].

I’d greatly appreciate any help you can give with this problem.

Many thanks,

Philip

PS: By the way, I can’t find any distinction in the code between ghost atoms that arise from parallelization and those that arise from periodic boundary conditions. Am I correct in thinking that both are treated the same?

Siemens Aktiengesellschaft: Vorsitzender des Aufsichtsrats: Gerhard Cromme; Vorstand: Peter Löscher, Vorsitzender; Johannes Feldmayer, Heinrich Hiesinger, Joe Kaeser, Rudi Lamprecht, Eduardo Montes, Jürgen Radomski, Erich R. Reinhardt, Hermann Requardt, Uriel J. Sharef, Klaus Wucherer; Sitz der Gesellschaft: Berlin und München; Registergericht: Berlin Charlottenburg, HRB 12300, München, HRB 6684; WEEE-Reg.-Nr. DE 23691322

2

Phillip,

All positions dispersions are run over comm.cpp. You will find an exchange etc. Looking at i.e. Verlet::integrate() will give you an idea of programmatic flow. Indeed, at every time step positions are communicated to neighboring ghost shells.

Currently none of the atom styles communicate velocities to the ghost shell. However, you can modify one of the atom styles to include velocities. You would have to include this in one of the atom_vec modules. I’d suggest copying atom_vec_full and create a new atom style. See style.h for integration and add to style_user.h

Pieter

3

Use the dpd atom style - those ghost atoms carry their velocity
with them, since it's used in the pair dpd calculations. Otherwise
it's the same as atom_style atomic.

Steve