How to connect molecule ID with compute results

LAMMPS LAMMPS Mailing List Mirror
1

Hey all,

I am trying to model simple chemical reaction by writing code where criteria for reaction will be 2 molecules reaching certain distance plus having enough energy to react . The problem is that fix bond/create only invokes reaction if 2 atoms reach certain distance, and need some additional code. It is OH anion which will make attack on C-atom , and energy of activation is aproximately kinetic+potential energy of OH anion. I can use compyre pe/atom and compute ke/atom to compute pe+ke energy of each anion in simulation , I can also calculate distance OH-Cattacked by compute distance/local and then use dump custom to dump : molecule ID , pe , ke , distance and write some external program to decide whether reaction will ocurr or not .

But what I dont know is if I use dump custom , will there be in same line : molecule ID , pe and ke(corresponding to that molecule ID) and distance ?

So my question is how can I connect molecule ID with results of compute ke/atom , ie to know which ke/atom result corresponds to which molecule ID , or if I called fix bond/create , how can I know which molecule ID paerticipated in bond creation ?

Anyone can help or give some advices ?

Thank you .

2

   Hey all,

  I am trying to model simple chemical reaction by writing code where
criteria for reaction will be 2 molecules reaching certain distance plus
having enough energy to react . The problem is that fix bond/create only
invokes reaction if 2 atoms reach certain distance, and need some additional
code. It is OH anion which will make attack on C-atom , and energy of
activation is aproximately kinetic+potential energy of OH anion. I can use
compyre pe/atom and compute ke/atom to compute pe+ke energy of each anion in
simulation , I can also calculate distance OH-Cattacked by compute
distance/local and then use dump custom to dump : molecule ID , pe , ke ,
distance and write some external program to decide whether reaction will
ocurr or not .

  But what I dont know is if I use dump custom , will there be in same line
: molecule ID , pe and ke(corresponding to that molecule ID) and distance ?

yes and no. you can output molecule id, pe/atom, ke/atom in a custom
dump but not distance/local, which would be only in a local style
dump.
even more so, since there is no compute distance/local.

   So my question is how can I connect molecule ID with results of compute
ke/atom , ie to know which ke/atom result corresponds to which molecule ID ,
or if I called fix bond/create , how can I know which molecule ID
paerticipated in bond creation ?

you can track the atom ids of atoms involved in explicit bonds via
compute property/local and you can track bonds with compute
bond/local. as local properties, these need to be output with dump
local. you can get atom id, molecule id, pe/atom, ke/atom all in a
custom style dump. it is up to you to write a program/script that
combines the information of those two dumps.

axel.

3

HI,

Thank you for your answer.

I can also compute distance by compute pair/local dist .

I have just one more question, after chemical reaction I will obviously have new molecule in simulation, so I need to define new atom types, pair_coeffs etc . Can I do it on beginning of simulation even if new molecule still doesnt exist or I do it right before bond is formed/broken?

Thank you.

4

All the pair coeffs are defined based on atom types, ditto

for bond coeffs based on bond types, etc. At the point

in the input script where you define them, it does not

matter if any actual atoms, bonds, etc exist with those

types.

Steve

5

HI,

Sorry I need to ask one more question to be sure about my program. As I said I need to monitor energies and decide which created bond to keep or reverse because fix bond/create enables to create bond only by distance criteria .

So my idea was following : I need to know when the new bond is created , one way I thought to do it is by compute property/local command :

group newbond type 4 14 #atom types which can create bond

compute 1 newbond property/local btype # c_1 will be zero until new bond is created

compute 2 newbownd property/local batom1 batom2 # ID of atoms which make new bond

So I thought I could check if first bond is created by :

if “c_1[1]==0” # or if “$a==0” whatever if correct by syntax

So my first question is if there is no bond is c_1[0] zero or NAN ?

My second question is this : lets say that 3 bonds are already created and accepted with energy criteria .

Next new bond I could check by : if “c_1[4]==0”

But can I be sure that latest bond will be c_1[4] ? Is it possible that newest bond will be placed in c_1[2] and old bonds will be in c_1[1] , c_[3] and c_1[4] ?

Thank you all .

6

You can’t do these kind of operations from a LAMMPS

input script. There are currently no commands that query the

status of one bond or let you break or create it.

There is also nothing on the “if” doc page that says it

allows syntax like you are writing, and even if

you used a variable, there is nothing on the

variable doc page that says a variable can have

access to a “local” compute like property/local.

Steve

7

Hi, thanx a lot for your answer . So now I know that I should definitely focuss on writting outter program and analyze data from dump files .