Ionization in LAMMPS -

LAMMPS LAMMPS Development
, ,
1

Hi there,

I’ve been using LAMMPS a lot for my research and it has been quite useful for simulating partially ionized plasmas !!

Now I’m interested in simulating a partially ionized plasma with variable ionization over time, this is, ionizing neutral atoms while running a simulation.

This is, the simulation starts with two species 1 and 2 using the same script I’ve always used (and it works) but this time I’ve modified the fix_nve.cpp in order to change the atom type from 2 to 1 (from neutral to ions), after the integration, with some hard-coded probability. The idea in the future is to include a MCC module with the ionization cross section if this simple test work.

I’ve only changed the final_integrate function in the fix_nve.cpp file as follows:

void FixNVE::final_integrate()
{
double dtfm;

// update v of atoms in group

double **v = atom->v;
double **f = atom->f;
double *q = atom->q;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;

if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}

} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];

    **if (type[i]==2){**

        **double r = dis(gen);**
        **double qe = 1.60217662e-19;**

        **if(r<0.00002){**

            **type[i] = 1;**
            **q[i] = qe;**


        }

    }

  }

}

atom->map_init();
atom->map_set();

}

and apparently it works… but not exactly as I’d like…

if I dump a file with the id, type and velocities of all the atoms I do see an decrease in #neutrals and an increase in #ions but at some point, after a few thousand timesteps I get this error:

Out of range atoms - cannot compute PPPM (…/pppm.cpp:1891)

Also if I compare the values of the kinetic energy computed by lammps for each species is different from the values I get computing it from the velocities for each species in the dump file.

I’m not sharing the script since is the same as before I did this change and it worked before, I’m probably not changing everything I need to in the .cpp file in order to “ionize” correctly or maybe I should write a separate fix file but I’m not sure on how to do it and what type of fix should it be…

If you can help me with this I’ll really appreciate it!

Thanks,
Marco

2

There are two concerns here:

  • Changing the charge of an atom can result in a significant change of force, which can make the simulation unstable and trigger errors like you have observed. You may thus need to change the charge in increments over multiple time steps to avoid issues.
  • Since you are using a long-range solver, you cannot just simply change some atom types and charges, you will also have to re-initialize the pre-computed structure factors of the long-range solver. Rather than using a (really ugly) hack for testing, as you are doing, you could just break your simulation into a series of short runs and then use LAMMPS’ scripting abilities to pick one or more suitable atoms and the set command to change its properties. With the next run command, LAMMPS will ensure that everything is properly initialized. Another long-range solver issue is that those expect the total system to be neutral. The more you are diverging from that assumption, the more questionable your results will be.
2 Likes
3

Hi! Thank you for your quick answer!

I understand, the problem is that I need to do a really long simulation and I need to have an ionization rate as continuous as possible. By doing it from the input script it will require hundreds if not thousands of lines of code to do this, that’s why I’m trying to avoid that (unless there is some simple way of doing it from the input script). If I turn-off the long range solver and increase the cutoff distance of the coulomb potential (within some acceptable value due to the computational cost), would that be enough?

Thanks!

4

You can do a loop. I am not suggesting this as a permanent solution. That would be creating a suitable new fix. But rather as an alternative to your integrator modification hack.

This is outside my area of expertise. You need to discuss this with somebody that has experience in the kind of simulations that you want to do. Or study the relevant published literature.

1 Like
5

Thank you for this tip!! I’ll try doing a loop, I think that can be useful until I figure out how to write an ionization fix correctly.

I’m sorry I think I didn’t formulate the question correctly. In the physics sense it’s okay as long as the pair correlation function is the same (which can be corroborated doing ocp simulations).

I mean if it is enough in the sense that only changing the charge and type of particles, lammps will recognize the new ions and compute the direct force (coul/cut) between them.

Thank you for your help!!

6

You can look through the code and see that (1) all Coulombic styles (including coul/cut) use the q[i] array to list charges (2) the commands you would use (such as set) change this array directly. So yes, charges you change using LAMMPS commands will be reflected in the pair and long-ranged charge forces.

Also be careful with units – in both real and metal units the electron charge is just 1.

Kolafa and Perram have an expression for the real-space truncation error in Ewald summation (https://www.tandfonline.com/doi/abs/10.1080/08927029208049126); you can use a similar approach to estimate the truncation error in using a pure 1/r^2 pair force. I have a feeling the “acceptable value” will turn out to be very imprecise.

2 Likes
7

In general truncating the Coulomb interaction (pair_style coul/cut) is bad even with a long cutoff, better would be a Wolf summation pair_style coul/wolf or pair_style coul/dsf

https://docs.lammps.org/pair_coul.html#pair-style-coul-dsf-command

1 Like