Hi,
I am trying to create a number rigid bodies on a lattice. I am using a molecule template defined by mySPH2.dat which is a simple dumbbell shape created with two spheroids with a 50% overlap. Using the following script I expect to get an assembly of 3375 (15x15x15) rigid molecules with 6750 atoms.
mySPH2.dat (160 Bytes)
Hi,
I am trying to create a number rigid bodies on a lattice. I am using a
molecule template defined by mySPH2.dat which is a simple dumbbell shape
created with two spheroids with a 50% overlap. Using the following script
I expect to get an assembly of 3375 (15x15x15) rigid molecules with 6750
atoms.################################
atom_style sphere
[...]
fix 1 all rigid/small molecule mol myMol
fix 2 all deform 1 x erate -0.01 y erate -0.01 z erate -0.01
run 1
##################################However running the code I get the following output:
LAMMPS (7 Dec
2015)
Read molecule
myMol:
2 atoms with 1
types
0 bonds with 0
types
0 angles with 0
types
0 dihedrals with 0
types
0 impropers with 0
types
Lattice spacing in x,y,z = 1 1
1
Created orthogonal box = (-11.25 -11.25 -11.25) to (11.25 11.25
11.25)
1 by 1 by 1 MPI processor
grid
*My ID: 0,
maxmol:0
*
*M**y ID: 0,
max:0
*
*M**y ID: 0,
molcreate:3375
*
*M**y ID: 0, moloffset:0 **Created 6750 atoms*
*1 rigid bodies with 6750 atoms*which suggests (bold underlined text) that only one body with 6750 atoms
is created. I put some additional print commands in create_atoms.cpp (lines
441-446) and seems initially 3375 molecules are created (italic text), I
have also attached a snapshot of the assembly, again 1 molecule is placed
within each of the 3375 unit cells. However, LAMMPS identifies them as a
single rigid object. Just wondering if this is the intended behavior, if so
how can I create 3375 rigid bodies each containing 2 atoms?
the first thing that comes to my mind is that atom style sphere does not
provide a molecule id, so fix rigid/small cannot detect molecules by
molecule id (as per "molecule" keyword). please note that the "mol" keyword
only is relevant for simulations where you are *depositing* molecules from
templates *during* the run (so that fix rigid/small knows about the rigid
body composition and properties).
thus try replacing "atom_style sphere" with "atom_style hybrid sphere bond"
and see if that helps.
axel.
And to verify you are creating molecules with create_atom, why
don’t you not use fix rigid and simply use a dump command with
the molID on step 0 and see what you have created.
Steve
Dear Steve and Axel,
Thank you for your comments. I tried axel’s suggestion (atom_style hybrid … ) but the code exits with an error (invalid style).
Also my impression from the documentation is that the “mol” keyword “can” be used to start a system with no rigid bodies initially (in conjunction with e.g. fix pour) but this is not a limitation?
The problem is create_atom correctly creates the molecules and sets the IDs (3375 molecules with 6750 atoms). So I think fix rigid should know about the molecule IDs (Mol Ids are added via fix property/atom mol and create_atom works correctly with this).
When I dump the atoms they all have a MolID=0 which is strange since the MolIDs get set in create_atom when they are created.
Many thanks
Sina
And to verify you are creating molecules with create_atom, why
don’t you not use fix rigid and simply use a dump command with
the molID on step 0 and see what you have created.
Steve
Dear Steve and Axel,
Thank you for your comments. I tried axel's suggestion (atom_style hybrid
... ) but the code exits with an error (invalid style).
Also my impression from the documentation is that the "mol" keyword "can"
be used to start a system with no rigid bodies initially (in conjunction
with e.g. fix pour) but this is not a limitation?The problem is create_atom correctly creates the molecules and sets the
IDs (3375 molecules with 6750 atoms). So I think fix rigid should know
about the molecule IDs (Mol Ids are added via fix property/atom mol and
create_atom works correctly with this).
When I dump the atoms they all have a MolID=0 which is strange since the
MolIDs get set in create_atom when they are created.
i've looked through the source code and i do *not* see where the molecule
id is set and thus getting molecule ids of 0 *does* make sense.
but in your case, there is a simple way to set them. after the
create_atoms command, simply add these two lines.
variable id atom ceil(id/2)
set atom * mol v_id
and then fix rigid/small should work as expected on the pairs of atoms. you
*still* should not use the mol keyword, because as far as i see, it cannot
compute the correct mass for rigid object of (overlapping) extended
particles. those would normally have to be read from a file. since you are
not depositing molecules during the run, the mol keyword is not needed.
axel.
Dear Axel,
Thank you for your suggestions, setting the molecule IDs manually did the trick and you’re right mol is not required in this case.
Just a quick question for me to better understand the code. I thought lines 434 onward in create_atom.cpp is an attempt to correctly calculate molecule ids for template style?
Many thanks
Sina
Dear Steve and Axel,
Thank you for your comments. I tried axel’s suggestion (atom_style hybrid … ) but the code exits with an error (invalid style).
Also my impression from the documentation is that the “mol” keyword “can” be used to start a system with no rigid bodies initially (in conjunction with e.g. fix pour) but this is not a limitation?The problem is create_atom correctly creates the molecules and sets the IDs (3375 molecules with 6750 atoms). So I think fix rigid should know about the molecule IDs (Mol Ids are added via fix property/atom mol and create_atom works correctly with this).
When I dump the atoms they all have a MolID=0 which is strange since the MolIDs get set in create_atom when they are created.
i’ve looked through the source code and i do not see where the molecule id is set and thus getting molecule ids of 0 does make sense.
but in your case, there is a simple way to set them. after the create_atoms command, simply add these two lines.
variable id atom ceil(id/2)
set atom * mol v_id
and then fix rigid/small should work as expected on the pairs of atoms. you still should not use the mol keyword, because as far as i see, it cannot compute the correct mass for rigid object of (overlapping) extended particles. those would normally have to be read from a file. since you are not depositing molecules during the run, the mol keyword is not needed.
axel.
Dear Axel,
Thank you for your suggestions, setting the molecule IDs manually did the
trick and you're right mol is not required in this case.
Just a quick question for me to better understand the code. I thought
lines 434 onward in create_atom.cpp is an attempt to correctly calculate
molecule ids for template style?
you are right, but there is a little problem with the logic. the code in
question will only apply those molecule id, in case you have a molecular
system with bonds and beyond. before it was possible to add the molecule
id via fix property/atom, this logic was sufficient. it should be possible
to update this section to be active also in this special case.
axel.
Dear Axel,
Thank you for your suggestions, setting the molecule IDs manually did the
trick and you're right mol is not required in this case.
Just a quick question for me to better understand the code. I thought
lines 434 onward in create_atom.cpp is an attempt to correctly calculate
molecule ids for template style?you are right, but there is a little problem with the logic. the code in
question will only apply those molecule id, in case you have a molecular
system with bonds and beyond. before it was possible to add the molecule
id via fix property/atom, this logic was sufficient. it should be possible
to update this section to be active also in this special case.
please try this change:
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index 951d1a7..19acd30 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -445,7 +445,7 @@ void CreateAtoms::command(int narg, char **arg)
tagint moloffset;
tagint *molecule = atom->molecule;
- if (molecule) {
+ if (atom->molecule_flag) {
tagint max = 0;
for (int i = 0; i < nlocal_previous; i++) max = MAX(max,molecule[i]);
tagint maxmol;
@@ -486,7 +486,7 @@ void CreateAtoms::command(int narg, char **arg)
for (int i = 0; i < molcreate; i++) {
if (tag) offset = tag[ilocal]-1;
for (int m = 0; m < natoms; m++) {
- if (molecular) molecule[ilocal] = moloffset + i+1;
+ if (atom->molecule_flag) molecule[ilocal] = moloffset + i+1;
if (molecular == 2) {
atom->molindex[ilocal] = 0;
atom->molatom[ilocal] = m;
axel.
Dear Axel,
Thank you, I will try this and let you know how it goes.
Many thanks
Sina
Dear Axel,
Thank you for your suggestions, setting the molecule IDs manually did the trick and you’re right mol is not required in this case.
Just a quick question for me to better understand the code. I thought lines 434 onward in create_atom.cpp is an attempt to correctly calculate molecule ids for template style?
you are right, but there is a little problem with the logic. the code in question will only apply those molecule id, in case you have a molecular system with bonds and beyond. before it was possible to add the molecule id via fix property/atom, this logic was sufficient. it should be possible to update this section to be active also in this special case.
please try this change:
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index 951d1a7…19acd30 100644
— a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -445,7 +445,7 @@ void CreateAtoms::command(int narg, char **arg)
tagint moloffset;
tagint *molecule = atom->molecule;
axel.
I’ll fix this in a patch today - good catch Axel.
Steve