Fix for particle creation from previous ones

Submitted by tkulju on Tue, 11/08/2011 - 15:24

Hi!
Now I'm trying to write the fix, which would create a new particle as particles are near to each other, i.e. A+B->C type of reaction. So would the best strategy be to take from fix_bond_create_gran.cpp the list building and implement it on the fix_insert_pack.cpp ? As I understood the atom positions are gone through in
void FixBondCreateGran::post_integrate()
and the insertion of the particles in
void FixInsertPack::x_v_omega
Am I correct? Or is there some more complicated things, which I should consider?

- Timo

Forums:

ckloss's picture

ckloss | Wed, 11/09/2011 - 09:52

Hi Timo,
x_v_omega(...) is used to determine position, velocity, omega of the particles to be inserted. So in your new FixInsertChemicalReaction : FixInsert you should loop a neighbor list in FixInsertChemicalReaction::x_v_omega(...) that is requested in FixInsertChemicalReaction::init(). Whenever a pair of particles is eligible for generating a third particle via a chemical reaction, increase ninserted_this and ninserted_spheres_this so that the mother class knows how many particles have been inserted

Cheers, Christoph

tkulju | Wed, 11/09/2011 - 15:12

Hi Christoph!
Thanks for the advice! So now when the particles react, how do I remove the old ones? I tried to read delete_atoms.cpp, but couldn't quite figure that out. Can I remove right after the suitable atoms are found in neighbor list loop, or do I have to search first all the particles, which could participate in the reaction?

I did the particle insertion with
ninserted_spheres_this+= pti->set_x_v_omega(pos,v_insert,omega_insert,quat_insert);. Does this do the job? For now it is unclear, how should I define v_insert,omega_insert,quat_insert for the particles to be inserted. Any suggestions?

Thanks,
Timo

ckloss's picture

ckloss | Thu, 11/10/2011 - 09:04

>>For now it is unclear, how should I define v_insert,omega_insert,quat_insert for the particles to be inserted. Any suggestions?
I am not a chemist, but I guess momentum should still be conserved?

>>So now when the particles react, how do I remove the old ones?
That's the bit trickier part to orchestrate all this so that it makes sense and is efficient. If you maybe have 2 weeks of patience, I would release a first version of a particle breakup model that would do something very similar (remove one particle and insert a couple smaller particles instead)

>>Can I remove right after the suitable atoms are found in neighbor list loop,
>>or do I have to search first all the particles, which could participate in the reaction?
I would first search through the neighlist, and if you find 2 particles eligible for reaction, store the properties of the new particle to insert in a list and delete the two particles that react. And then, after looping the neighlist, loop this list and generate all the new particles.

Cheers, Christoph

tkulju | Thu, 11/10/2011 - 09:43

Hi!
Ok, I think the momentum conservation would be a good criteria. Can you give any hints, how I can remove the particles? Or is this a more complicated issue?

Thanks,
Timo

ckloss's picture

ckloss | Thu, 11/10/2011 - 11:03

>>Can you give any hints, how I can remove the particles?
If you can wait for the next release next week, you will see how it is done in combination with particle insertion in fix breakparticle/force

Cheers, Christoph

tkulju | Thu, 11/17/2011 - 15:14

Hi!
Now I'm getting some understanding how this fix insert stuff works. When I do some modification to my fix, do I have to compile the whole code again (via make clean-all; make fedora) ? Or is there any quicker way to get the code implemeted in the already existing lmp_fedora?

- Timo

dbreton's picture

dbreton | Mon, 11/21/2011 - 17:49

Yes, you need to rebuild all of LIGGGHTS every time you want to implement a change you made in the source code. The first re-build will take some time, but subsequent ones are faster, assuming you are only changing small bits of code. The -j switch in make can speed things up a little by using more than one processor (i.e. make -j4 fedora), but I find that it does not make much difference... check out the man page if interested.

good luck
DAN

tkulju | Fri, 11/18/2011 - 09:57

Hi!
Ok, how can I tell the program that "don't try to insert any particles, if ninserted_this==0"? I guess I have to do this in FixInsertReaction::init(), after FixInsertReaction::x_v_omega(...) has been called, but don't know how to pass this message to LIGGGHTS...

Otherwise the program will stop with an error "ERROR: Illegal fix insert command, would insert <=1 particles per insertion step", if the ninsert_per<=1.

- Timo

tkulju | Mon, 11/21/2011 - 07:53

Figured that out; plain return-command does the job.. I guess inserting it directly to x_v_omega(...) would be the wisest thing to do?

- Timo

tkulju | Mon, 11/21/2011 - 10:57

Hi!
Now I have a problem with the lists. When I'm tying to get the firstneighbor (?) via list->firstneigh, I'm getting a segfault. But for other list variables, such as list->inum, list->ilist, list->numneigh I can get just fine. Any idea, what could be the problem? Do I have to allocate this firstneigh list first somehow?

- Timo

tkulju | Mon, 11/21/2011 - 11:54

My bad, the problem wasn't in the list but in some 2D pointer allocation. Found it and trying to fix it now...

- Timo

ckloss's picture

ckloss | Mon, 11/21/2011 - 12:55

Hi Timo,

there will be a slightly revised fix insert in the new release this week which makes it possible to define how many particles in the derived class. That would make it easier for you...so stay tuned!

Cheers, Christoph

tkulju | Tue, 11/22/2011 - 09:24

Hi!
Before I can use the lists, I think I have to allocate them first? If so, how can I do that? I came to this conclusion by getting several segfaults, and finally adding a line in x_v_omega(...) (before the list calls)
if(list==NULL)error->all("...");
stopped the execution. Or is there something else, why I'm getting this error?

- Timo

ckloss's picture

ckloss | Tue, 11/22/2011 - 10:03

Did you request a list as done in compute pair/gran/local? then the list should not be null

Christoph

tkulju | Tue, 11/22/2011 - 12:03

Hi!
Ok, what I understood, I have to add neighbor requests to init() for this fix? When adding the following:

int irequest = neighbor->request((void *) this);
// if(screen)fprintf(screen,"irequest=%d\n",irequest);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;

if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;

I'm getting the segfault right after "Setting up run..." message. Otherwise it continues to the x_v_omega(...) (first called in pre_insert() in the beginning of FixInsert::pre_exchange()). I have the init_list(...) as in compute_pair_gran_local.cpp, which should initialize the list, am I right?

- Timo

ckloss's picture

ckloss | Tue, 11/22/2011 - 12:39

That's now quite a lot of technical details... I suggest you first accommodate yourself with the way a fix or compute can use a neighbor list, see eg fix bond/create or compute pair/gran/local

Cheers, Christoph

tkulju | Fri, 11/25/2011 - 12:31

Hi Christoph!
What I've understood is that in order to obtain the neighborhood lists, first I have declare the variable *list in header file as class NeighList *list. Then I can use it to get the needed variables from NeighList class. But I get list==NULL, if I don't declare the list to point to NeighList class in x_v_mega(...). And then I can't get the ilist (it crashes, when I'm trying to use it). I tried to initialize it as list->grow(nmax), but it doesn't help. So what am I missing? I couldn't find any further clues from fix_bond_create_gran.cpp or compute_pair_gran_local.cpp files, which you mentioned.

- Timo

ckloss's picture

ckloss | Fri, 11/25/2011 - 12:40

Timo,

you have to request a neighbor list (be sure to request a granular list with fix as requestor) and then the value is assigned to the class data member by means of a callback, have a look at
FixBondCreateGran::init()
and
FixBondCreateGran::init_list()

If this does not help, one more clue: As x_v_omega(... is called of pre_exchange(), I think but am not 100% sure if the neigh list is still valid there - I never tried that. You should try out at which point of the verlet time stepping (see developers manual) the list is available and then process it at the right point.

Cheers, Christoph

tkulju | Fri, 11/25/2011 - 14:22

Hi Christoph!
Thanks for the advise. Now the code progressed to the init_list phase, but after that it crashes. What I've done, I have set the neighbor->requests[irequest]>fix=1 (and the respective pair=0), where irequest=request((void *) this), as it was also in FixBondCreateGran::init().

The crash actually happens in pre_exchange at FixContactHistory, and I think it's a caused consequence when i=ilist[ii] (at fix_contact_history.cpp, line 176) gets a huge negative value. Any idea how to fix this?

- Timo

ckloss's picture

ckloss | Mon, 11/28/2011 - 12:38

hi, without knowing what exactly you have coded, I can give you two more hints:

+ try to use a memory debugger such as valgrind (compiling with -D will enable to to get the code line that causes the trouble)
+ Instead of requesting a neigh list on it's own, you can also access a neigh list that exists for sure - this is e.g. done in fix heat/gran, which accesses the pair gran list directly. Maybe not a very elegant way, but if you are not 100% sure what is the problem, you can try it that way

Hope I could help - cheers, Christoph

tkulju | Tue, 11/29/2011 - 14:30

Hi Christoph!
Thanks for your help, the latter one worked!! I didn't even try the first one... Now I can access the neighbor lists and the code doesn't crash!

- Timo

tkulju | Wed, 11/30/2011 - 15:27

Hi!
Now I have succesfully been able to run my fix for the first timestep. When the fix is executed at second timestep, the list variables do not seem to be updated. When deleting the particle, the last atom of the atom vector is inserted in the position of the deleted atom and the atom->nlocal is updated. Nothing is done for the neighbor (list) variables. Should I udate the neighbor (lists) also and if so, how?

Thanks,
Timo