LAMMPS library: Gather specific atoms and insert small molecules

    3) create_atoms can create/distribute
    random atoms. I need however to create
    small molecules and distribute them
    randomly. The small molecule structure
    is known a-priori, however I need to
    insert let's say 10 randomly in my
    simulation box.

with a recent enough version of LAMMPS, ​create_atoms can also create
molecules that have been previously defined with the molecule command
and read from a suitable file.

Follow up on this:
If I would add a new ion type to my system,
I would have to specify obviously all the
interactions with the other types in my system.
Correct? So would it be better to prepare the
ion in a let's say dummy PDB file and convert
it to a LAMMPS input data with e.g. topotools
and then use read_data?

Best regards,
Stephan

>
> 3) create_atoms can create/distribute
> random atoms. I need however to create
> small molecules and distribute them
> randomly. The small molecule structure
> is known a-priori, however I need to
> insert let's say 10 randomly in my
> simulation box.
>
>
> with a recent enough version of LAMMPS, ​create_atoms can also create
> molecules that have been previously defined with the molecule command
> and read from a suitable file.

Follow up on this:
If I would add a new ion type to my system,
I would have to specify obviously all the
interactions with the other types in my system.
Correct? So would it be better to prepare the
ion in a let's say dummy PDB file and convert
it to a LAMMPS input data with e.g. topotools
and then use read_data?

​the read_data command takes options like: ​extra/atom/types
where you can reserve space for these atom types.

the basic rules in LAMMPS related to this are simple:
- number of atom/bond/angle/dihedral/improper types per atom, max
bonds/angles/dihedrals/impropers per atom, max special per atom are locked
in, when the simulation cell is created. this happens with either
create_box, read_data or read_restart.
- this cannot be changed after that as it would require to redo the box
creation. a new box can only be created after the clear command, which
reinstantiates the LAMMPS class from scratch.

​no need to use a dummy atom, tho.​

​axel.​

Is there a way to get the first extra atom type ID?
(Or: Get the numbers of used atom type IDs?)

Then I could specify pair coeffs etc. for the
newly added atom types.

Best regards,
Stephan

Is there a way to get the first extra atom type ID?
(Or: Get the numbers of used atom type IDs?)

a) ​you know how many extra types you requested.
b) you can query the total number of types present after the box has been
created.​
a) + b) => ntypes - nextra + 1 = first new atom type id.

<insert snarky comment here>,
     axel.

    Is there a way to get the first extra atom type ID?
    (Or: Get the numbers of used atom type IDs?)

a) ​you know how many extra types you requested.
b) you can query the total number of types present after the box has
been created.​
a) + b) => ntypes - nextra + 1 = first new atom type id.

Yes, this is true.
Thank you.
(Using lmp->atom->ntypes should be fine).

Best wishes,
Stephan

    >
    > 3) create_atoms can create/distribute
    > random atoms. I need however to create
    > small molecules and distribute them
    > randomly. The small molecule structure
    > is known a-priori, however I need to
    > insert let's say 10 randomly in my
    > simulation box.
    >
    >
    > with a recent enough version of LAMMPS, ​create_atoms can also create
    > molecules that have been previously defined with the molecule command
    > and read from a suitable file.

    Follow up on this:
    If I would add a new ion type to my system,
    I would have to specify obviously all the
    interactions with the other types in my system.
    Correct? So would it be better to prepare the
    ion in a let's say dummy PDB file and convert
    it to a LAMMPS input data with e.g. topotools
    and then use read_data?

​the read_data command takes options like: ​extra/atom/types
where you can reserve space for these atom types.

the basic rules in LAMMPS related to this are simple:
- number of atom/bond/angle/dihedral/improper types per atom, max
bonds/angles/dihedrals/impropers per atom, max special per atom are
locked in, when the simulation cell is created. this happens with either
create_box, read_data or read_restart.
- this cannot be changed after that as it would require to redo the box
creation. a new box can only be created after the clear command, which
reinstantiates the LAMMPS class from scratch.

​no need to use a dummy atom, tho.​

Thanks for pointing that out.

This might seem like a silly question:
If I don't need extra atom types for my molecule
I want to add via the molecule command I could just
use the atom types provided from my "main molecule"
which I loaded before with the read_data command.
There is probably no automation I could make use of.
Correct? (Just want to make sure I don't miss
something or do unneccesary work). TIA.

Best wishes,
Stephan

>
>
>
> >
> > 3) create_atoms can create/distribute
> > random atoms. I need however to create
> > small molecules and distribute them
> > randomly. The small molecule structure
> > is known a-priori, however I need to
> > insert let's say 10 randomly in my
> > simulation box.
> >
> >
> > with a recent enough version of LAMMPS, ​create_atoms can also
create
> > molecules that have been previously defined with the molecule
command
> > and read from a suitable file.
>
> Follow up on this:
> If I would add a new ion type to my system,
> I would have to specify obviously all the
> interactions with the other types in my system.
> Correct? So would it be better to prepare the
> ion in a let's say dummy PDB file and convert
> it to a LAMMPS input data with e.g. topotools
> and then use read_data?
>
>
> ​the read_data command takes options like: ​extra/atom/types
> where you can reserve space for these atom types.
>
> the basic rules in LAMMPS related to this are simple:
> - number of atom/bond/angle/dihedral/improper types per atom, max
> bonds/angles/dihedrals/impropers per atom, max special per atom are
> locked in, when the simulation cell is created. this happens with either
> create_box, read_data or read_restart.
> - this cannot be changed after that as it would require to redo the box
> creation. a new box can only be created after the clear command, which
> reinstantiates the LAMMPS class from scratch.
>
> ​no need to use a dummy atom, tho.​

Thanks for pointing that out.

This might seem like a silly question:
If I don't need extra atom types for my molecule
I want to add via the molecule command I could just
use the atom types provided from my "main molecule"
which I loaded before with the read_data command.
There is probably no automation I could make use of.
Correct? (Just want to make sure I don't miss
something or do unneccesary work). TIA.

​correct. there is no way to automate this. by the time you could check
whether there are redundant atom types through comparing pair_style and
pair_coeff settings, the number of atom types is already locked in. for
most simple force fields, you only need to provide the per type parameters
and then the mixed parameters can be generated from the selected mixing
rule. anything beyond that (e.g. many-body or hybrid pair style setups, you
will have to figure it out "manually" on a case by case basis.

axel.

correction. on second thought, i think that you could do some processing of settings and parameters.
in that case you would have to set up your system once, then identify redundant settings and store this information, delete everything with “clear” and start all over using the information from the previous setup attempt.

but in my opinion, it is not worth the hassle. i would just set up everything assuming no overlap of types and then script that process. it seems much easier and simpler to me.

axel.

I think I will go with that...

but in my opinion, it is not worth the hassle. i would just set up
everything assuming no overlap of types and then script that process. it
seems much easier and simpler to me.

I managed to set up my force field, though, now I
inserted by create_atoms my molecules randomly.

Is it correct to use (if my molecule has ID 1)
"group MYGRP molecule 1" to get all the replicated
molecules in my box?

Because then I could use the code you gave me
in a previous answer to extract all coordinates
from these molecules by:

lammps_command(lmp,"variable MYCOORDSX atom x");
                                                                    double *coordsX = (double
*)lammps_extract_variable(lmp, "MYCOORDSX","MYGRP");

This seems not to give me all x-coordinates
(I replicated 3 molecules with 19 atoms but
I'm getting only 6 x-coordinates.)

Thanks in advance.
                    
Best regards,
Stephan

Dear LAMMPS users,

let me rephrase my previous post to be more precise.
I read in a protein structure by and added
after that three molecules randomly by:
read_data main_protein.data
molecule 1 small_protein.data
create_atoms 1 random 3 1506024437 NULL mol 1 91731

Then, I want to retrieve the atoms of the three inserted
molecules (which have a total of 19 atoms each) by:
group MYGRP molecule <> 1 3

This will however not give me 57 atoms (3*19 atoms):
48 atoms in group MYGRP

I think I'm missing a piece of information
(Or more likely doing something wrong...)

Thanks in advance.

Best regards,
Stephan

I suggest you dump a snapshot and see what

the molecule IDs are of the atoms you added.

Steve

Hey Steve,
thanks for your suggestion.

This revealed that my molecules are created
with molecule ids 9260, 9261, 9262.
Now the question would be: If I know how
many molecules I'm creating, how would I
get the first index, e.g. 9260 by the
LAMMPS library?

Best regards,
Stephan

Hey Steve,

I'm defining my group as follows:
Determine the first index (by looking at the dump) of
my molecule of interest and calculate last index by:
last index = first index + nMols.
(nMols I know because I added e.g. nMols = 3 molecules).

Then I can form a group.
I found out that this question has been asked
on the mailing list (archive) but I could not
find an answer. That's why I politely ask again:

How can I get the total number of molecules
in my simulation box (Preferably by a call to
lammps_command)? Thanks in advance.

Best wishes,
Stephan

Addendum: I would also be willing to implement
a library function for that purposes.

Best wishes,
Stephan

There is no LAMMPS command that will tell

you the current # of molecules. Since molecule IDs

do not need to be contiguous, I can’t think of a simple

way to calculate it. E.g. by a compute.

You could use compute reduce to calculate the max or

min molecule ID. Which means you could assign the

result to a variable and invoke the variable via

the lib interface, e.g. lammps_command() or

lammps_extract_variable(). If you figure out

how to count unique molecule IDs (in parallel), you could

write it yourself external to LAMMPS and call

extract_atom() thru the lib interface to get pointer

(on each proc) to the mol IDs.

Steve