Compute displace/atom & electric ionic current

LAMMPS LAMMPS Beginners
1

Dear Lammps Users,
I want to calculate instantaneous current, as mentioned below.

Screenshot 2022-07-25 162747
(https://pubs.acs.org/doi/10.1021/acsanm.9b02280?ref=pdf)

Initially I thought of using compute displace/atom command and as a result I went through its documentation.

Define a computation that calculates the current displacement of each atom in the group from its original (reference) coordinates, including all effects due to atoms passing through periodic boundaries.

A vector of four quantities per atom is calculated by this compute. The first 3 elements of the vector are the dx,dy,dz displacements. The fourth component is the total displacement, i.e. sqrt(dxdx + dydy + dz*dz).

The displacement of an atom is from its original position at the time the compute command was issued. The value of the displacement will be 0.0 for atoms not in the specified compute group.

My rationale was since a_i(t + \Delta t) - a_i(t) is displacement hence I could use the command here. At the moment, I not only don’t know how to calculate instantaneous current but also I’ve trouble understanding how compute displace/atom works.

Suppose I have created a group with lets say 10 atoms all of same type:

group         abc id 1:10
compute       disp abc displace/atom 
thermo_style  custom step time c_disp

I understand this will give me an error because because what I’m computing is a n * 4 dimensional vector and what I’m trying to print is a scalar but even so my question is how this command will print out displacement of each atom in the group? In other words, how can I write it in a way thru which I get displacement of each atom along a particular direction, say x-axis? More importantly, can I make use of this command for my problem?

Even if I use something like

group         abc id 1:10
compute       disp abc displace/atom 
compute       disp2 abc reduce ave c_disp[1]
thermo_style  custom step time c_disp

Now disp, which is a vector, has been converted to scalar using reduce command. But it’d be blasphemous to use this to calculate the above mentioned instantaneous current problem.

I don’t know how to make use of trajectory file generated after a particular LAMMPS simulation for such a calculation. Any help regarding how can I make use of trajectories in calculating the above mentioned current problem would be much appreciated.

Best
Akshay

2

First of all. You are mixing up total displacement, i.e. a(t) - a(0), with per step displacement, i.e. a(t+delta t) - a(t)).

Second, the formula you quote is the sum of the individual components multiplied by their charge. So using the fourth component or the average in a reduction is clearly not consistent with that.

It is very dangerous to just try out commands and hope that they produce the properties you want based on just speculation and the fact that LAMMPS does not reject it. You have to figure out the math and for that you have to study the corresponding documentation (and if that is not sufficient) the corresponding source code.

If you want the accurate displacement between two timesteps you could use fix store/state to store the (unwrapped!) positions of the previous step, then define an atom style variable to compute the displacement times the charge and then use compute reduce sum to get the total and then multiply that result as indicated in the equation you quote with the suitable prefactor.

1 Like
3

You could use the refresh option to reset the origins every few time steps, but I have no personal experience with that. In your situation I would do:

group       typeA        id 1:10 # descriptive names help you remember stuff!
compute     atomdispA    typeA   displace/atom
compute     groupdispA   typeA   compute/reduce sum atomdispA[*]
fix         outputdispA  all     print 100 &
       "$(c_groupdispA[1]) $(c_groupdispA[2]) $(c_groupdispA[3])" & 
       file outputdispA.txt

Note that the compute/reduce here reduces a 4-column per-atom variable to a 4-value global vector. In post-processing, the values you are looking for will be the differences between successive rows, multiplied by a suitable factor. Note that this assumes all particles in the group typeA have the same partial charge.

Most importantly, note that you have to make sure you understand all this for yourself! If you are asked what this code does (including by future you), “a random stranger told me this on the Internet” is a deeply unsatisfying explanation.

1 Like
4

If I look at the quoted expression, it looks like it is trying to reconstruct the velocity from delta x by delta t. But in LAMMPS you can access vx directly and thus the current in x would be the sum of q*vx divided by the width of the slice you are summing over. Intuitively this would make sense to me, but this is not my area of research.

I would look for a proper derivation of the property you are after in a credible text book.

1 Like
5

From what I know, if an average flow is all that’s needed, it’s more statistically efficient to use displacements. The displacement every N time steps automatically integrates N “samples” of the velocity and thus the variance is reduced (a little, since those samples are correlated, but more and more as N increases). See https://pubs.acs.org/doi/10.1021/jp5012523 for more advanced treatments.

Of course, if correlations or fluctuations in velocity are important, then the actual velocities do need to be sampled (although these expressions can often be analytically integrated into displacement measurement versions). Also, in periodic simulations with changing box lengths, particle displacements and velocities are a combination of “regular motion” and box length changes, which may need to be separately treated – this affects measurements of diffusivities from MSDs, see Systematic errors in diffusion coefficients from long-time molecular dynamics simulations at constant pressure: The Journal of Chemical Physics: Vol 153, No 2.

1 Like
6

What I understood is that its the sum of any one coordinate (that is decided by the direction in which electric field is applied; suppose x axis so, that’d be the sum of x_i(t+\Delta t) - x_i(t)) for all the N atoms multiplied by their individual q's.

I’m a bit confused; can you please elaborate the difference between the two a bit more?

I’ll look for the derivation as well. Thanks for the suggestion.

@srtee Thanks for these two papers. I’ll go thru them and try to figure out how can the concepts be applied to my problem.

7

Dear Lammps User,
My system consists of 300 ions in total. 150K+ and 150 Cl-. A polymer chain of 50 beads whose backbone is modeled by FENE potential. And lastly, there is membrane with a pore in it. I tried to make changes as suggested by @akohlmey & @srtee but due to my naivete I’m not being able to implement. The primary idea is still the same which is to calculate

I tried and saw that the displace/atom command is not producing the results which I should have gotten according to the theory which I’m trying to model. As a result, I am going for fix store/state. The rationale is to fix the state (i.e. position) of each bead for let’s say t time and then again fix it for t+ \Delta t and further by subtracting the two and summing (compute/reduce) it over all ions, should yield me the correct results.

However, I’m not sure how to go about it?

I tried writing the below code but I think I’m butchering the logic and not going anywhere with it.

fix prev_pos ions store/state 2 xu yu zu
fix curr_pos ions store/state 1 xu yu zu
variable pos_x equal “(f_prev_pos-f_curr_pos)”

and then further multiplying with the mentioned prefactors and then dumping out these values.

Can you please help me figure out whether I’m headed in the right direction or not?

Best
Akshay

8

Why not work from the dump files and use post-processing to calculate the quantity you’re looking for?

1 Like
9

I don’t think that there is any functionality currently in LAMMPS that can be used in the way you intent to use it. What you would need is the per-atom equivalent to the fix vector command. This has currently been avoided because of the excessive memory needs if you do not limit the maximum vector length.

However, I have to come back to my statement from earlier: Since the quoted segment explicitly is stating that it is computing the displacement for a single(!) timestep, I maintain that using the velocity vx instead of (x(t+delta t) - x(t))/delta t is the best approach.

variable lx index 10.0
variable ly index 10.0
variable lz index 10.0
variable ex atom q*vx
variable ey atom q*vy
variable ez atom q*vz

compute current all compute reduce sum v_ex v_ey v_ez
variable cx equal c_current[1]/v_lx
variable cy equal c_current[2]/v_ly
variable cz equal c_current[3]/v_lz

cx, cy, cz would then be the current in the corresponding directions (there is no mention of using the total displacement), these can then be fed to a fix ave/time command to get less noisy data averaged over multiple steps.

1 Like
10

I’m running these for 8e+6 iterations with 0.005 t^* timesteps, equivalent to 18.5 \mu s. And this I have to repeat 250-300 times for averaging. The post processed file is way to big and it takes a bit longer to plot the data on my desktop. Whereas I can run the same code using loop command in LAMMPS using my institute’s cluster computers and dump out only the relevant information in a separate file and use that data for further analysis. This is the way I have been doing up until now, If there is a room to fasten this process I’ll discuss it with HPC people in my institute.

I didn’t look for the derivation earlier, using a textbook but I’ll get on with it now. Thank you so much @akohlmey .