Simulation of Molecule Aggregation in a Vesicle Using LAMMPS

Simulation of Molecule Aggregation in a Vesicle Using LAMMPS

Hello everyone,

I would like to simulate the aggregation process of several coarse-grained spheres in a vesicle using LAMMPS.

The tests I plan to conduct are as follows:

  1. Vesicle only.
  2. Vesicle + water (represented by one sphere).
  3. Vesicle + water + negatively charged coarse-grained spheres. I hope to see no aggregation.
  4. Vesicle + water + negatively charged coarse-grained spheres + positively charged coarse-grained spheres. I hope to see aggregation.

For testing purposes, I will only use 2 coarse-grained atoms of each type for now.

  • Test 01 (Vesicle only): For this case, I referred to an example in LAMMPS under examples/ASPHERE/vesicle and was able to complete the simulation normally.
  • Test 02 (Vesicle + Water): I added water spheres to the system (vesicle with 2 atoms and two water spheres for testing). However, I encountered an error: “atom format error. All pair coeffs are not set (…/pair.cpp:253).”

Since I am using the YLZ potential, the vesicle needs to be defined as an ellipsoid.

My questions are:

** 1. For the negatively charged spheres, which atom type(s) should I use? (I am currently using charge.)
** 2. I have researched the solution to the “All pair coeffs are not set” error. It appears to be caused by the following pair definitions:

pair_style hybrid ylz 2.6 lj/cut 3.6
pair_coeff 1 1 ylz 1.0 1.0 4 3 0.0 2.6
pair_coeff 1 2 lj/cut ${ew} 1.0
pair_coeff 2 2 lj/cut ${ew} 1.0

To resolve this, I tried the following:

pair_style hybrid ylz 2.6 lj/cut 3.6 lj/cut 3.6
pair_coeff 1 1 ylz 1.0 1.0 4 3 0.0 2.6
pair_coeff 1 2 lj/cut 1 ${ew} 1.0
pair_coeff 2 2 lj/cut 2 ${ew} 1.0

However, this didn’t work. Can anyone provide some guidance on how to fix this?

Thank you

test.in

units lj

variable T equal 0.23
variable     ew  equal 0.2

atom_style hybrid ellipsoid charge

boundary     p p p  

read_data read_data.atoms

mass * 1.0

group mem type 1
group wat_in type 2

velocity mem create ${T}  87287 loop geom 
velocity wat_in  create ${T}  87287 loop geom 

pair_style   hybrid ylz 2.6 lj/cut 3.6 
pair_coeff   1  1 ylz 1.0  1.0  4  3  0.0  2.6
pair_coeff   1  2 lj/cut  ${ew} 1.0 
pair_coeff   2  2 lj/cut  ${ew} 1.0 

neighbor    1.0 bin

thermo_style custom step temp press ebond 

thermo 200
timestep 0.005

fix 1 wat_in npt temp 0.23 0.23 1 iso 0.05 0.05 1
run 100


read_data.atoms

vesicle with diameter =    30.000000
       4 atoms
       2 ellipsoids

5 atom types

  -35.000000    35.000000 xlo xhi
  -35.000000    35.000000 ylo yhi
  -35.000000    35.000000 zlo zhi

# type 1 : membrane CG,    ellipsoids
# type 2 : Water in Cell,  atomic
# type 3 : Water out Cell, atomic
# type 4 : dna,            atomic + charge
# type 5 : prob,           atomic + charge

Atoms  #atom ID, type,     x,            y,           z,         ell flg,   density, charge              
 
       1      1            0.236369      0.008147     0.007400   1         0.238732   0.0 
       2      1            0.236369      0.009058     0.008176   1         0.238732   0.0
       3      2            1.236369      0.008147     1.007400   0         0.238732   0.0
       4      2            0.236369      1.008147     1.007400   0         0.238732   0.0
#       5      3            2.236369      0.008147     1.007400   0         0.238732   0.0
#       6      3            0.236369      2.008147     1.007400   0         0.238732   0.0
#       7      4            3.236369      0.008147     1.007400   0         0.238732   -2.0
#       8      4            0.236369      3.008147     1.007400   0         0.238732   -2.0
#       9      5            4.236369      0.008147     1.007400   0         0.238732   2.0
#      10      5            0.236369      4.008147     1.007400   0         0.238732   2.0


Ellipsoids

       1     1.000000     1.000000     1.010000      0.707299     0.000000    -0.706915     0.000349
       2     1.000000     1.000000     1.010000      0.707320    -0.000000     0.706893     0.000385


By the way. I don’t know why the background color of test.in is black.So weired.

You are talking about the atom style and not the atom type here. The style “charge” will work, but it will be rather pointless to assign charges, when neither or your pair styles actually uses that information to compute forces.

Your data file lists 5 atom types, but you provide pair_coeff settings covering only 2 of them. It is no surprise that LAMMPS is telling you that there are settings missing. You can easily cross-check this by adding an “info coeffs” command.

Not for me.

@akohlmey

Thank you so much for the replies. They help me a lot.

You are talking about the atom style and not the atom type here.

Thank you so much. My BAD.

The style “charge” will work, but it will be rather pointless to assign charges, when neither or your pair styles actually uses that information to compute forces.

Since I will change the size of the CG in the future. I think sphere + charge may be the right option. Am I right?

Your data file lists 5 atom types, but you provide pair_coeff settings covering only 2 of them. It is no surprise that LAMMPS is telling you that there are settings missing. You can easily cross-check this by adding an “info coeffs” command.

So, that is the reason. And the “info coeffs” works.

It is not possible to answer this question in this generality. It all depends on the pair style. When your pair style is lj/cut, for example, neither atom style sphere nor atom style charge are more useful than the default of atom style atomic, since the lj/cut pair style will neither use the radius information nor the charge information. For the Lennard-Jones potential, the effective radius of the particles is determined by the sigma parameter of the pair coefficient i.e. the potential directly.

In short, it all depends on the model you are simulating and the potential you are using for that.

1 Like

@akohlmey
First, thank you so much for you kindness reply.

According to your suggestion, I changed LJ to LJ/cut/coul/long to take into account the electrostatic interaction between charges.
To verify this setting, I build a simple example. I wanted to simulate the interaction between an ellipsoid with a positive charge and an ellipsoid with a negative charge, in order to observe the attraction between the two coarse-grained particles. The input file is as follows.

The results show that the two particles indeed attract each other and cluster together, but then they quickly separate and continue to oscillate between attraction and repulsion. I am not sure where the issue lies. From my understanding, once the two ellpisoids aggregated, they should not separate again, but instead move as a single unit.

I plotted the distance between the centers of the two ellipsoids, and at the closest point, the distance between the centers is only 1.2. However, each ellipsoid should have its own radius. The distance between the centers of the two ellipsoids should be greater than the sum of their radii. This is also what confuses me.

Thank you for your reply.

# Initialize simulation
units lj
dimension 3
boundary p p p
atom_style hybrid ellipsoid charge

# Create simulation box
region box block -10 10 -10 10 -10 10
create_box 2 box

# Create positively charged sphere and negatively charged ellipsoid
create_atoms 1 single 0 0 0
create_atoms 2 single 5 0 0

# Set particle types
set atom 1 type 1  # Sphere
set atom 2 type 2  # Ellipsoid

# Set sphere properties
set atom 1 shape 1.0 1.0 1.1
set atom 1 charge 1.0
set atom 1 mass 4.18879  # Mass = (4/3)*pi*r^3 * density (assuming density is 1)
set atom 1 quat 1.0 0.0 0.0 0.0  # Initial orientation (quaternion representation)

# Set ellipsoid properties
set atom 2 shape 2.0 1.0 1.0  # Ellipsoid shape (a=2.0, b=1.0, c=1.0)
set atom 2 charge -1.0
set atom 2 mass 8.37758  # Mass = (4/3)*pi*a*b*c * density (assuming density is 1)
set atom 2 quat 1.0 0.0 0.0 0.0  # Initial orientation (quaternion representation)

# Define potential functions
pair_style lj/cut/coul/long 10.0
pair_coeff 1 1  1.0 1.0 2.5  # Lennard-Jones interaction between sllipsoids 1
pair_coeff 2 2  1.0 1.0 2.5  # Lennard-Jones interaction between ellipsoids 2
pair_coeff 1 2  1.0 1.0 2.5    # Coulomb interaction between ellipsoid 1 and ellipsoid 2

# Set long-range Coulomb interaction
kspace_style pppm 1.0e-5

# Set initial velocities
velocity all create 1.0 4928459 rot yes dist gaussian

# Run simulation
fix 1 all nve
thermo 100
thermo_style custom step temp pe ke etotal
dump 1 all atom 100 dump.traj.lammpstrj
run 10000

Here is the issue. You seem to lack a sufficient understanding of elementary classical mechanics.
If the atoms would just stick together, where would their kinetic energy go?

You have no ellipsoids in this simulation. Your pair style does not pay any attention to the ellipsoid settings.

Overall, you seem to be insufficiently prepared to do what you do on your own. Your questions strongly hint that you need some training by a person that has sufficient knowledge. Doing this on your own is going to lead to more problems and numerous mistakes, since you are already strugging at an extremely basic level. An online forum is no substitute for training. The people here cannot fill in for being your tutor. We can point you in the right direction, but the rest is your task and if you not (yet) up to it, you must get sufficiently competent local assistance.

You already see in our conversation, that despite my previous explanations, you have made the same mistake as before, only reversed, because you have not fully understood what I told you.

@akohlmey Thank you so much for your kind suggestions.
Since there is no local people who can advise me or discuss with me, I have to post my question and confusion. I will get the basic level of knowledge as you suggest.
Still, really thank you so much for these suggestions; though they did not solve my problem, but they point me to the “Right Direction”.

It is not a smart choice to start a research project in a new field without having proper tutoring and supervision. You are bound to make many unnecessary mistakes and waste a lot of time and then there is a significant risk in that your results may be bogus without anybody noticing.

Besides, there are more people around with the knowledge you need, you probably have not yet looked in the right places or may need to reach out (or your adviser or supervisor) to start a proper collaboration. The fact that computers and simulation software are easily available and the software has matured to a point that one can run it without much knowledge, does not translate that leads to successful results. If you were to use an expensive instrument (say a new and expensive spectrometer), nobody would even consider letting you start a project without proper training. Yet, with computer simulations people think they can manage because there is no cost involved.

Online forums are no replacement for proper training. There is no problem is asking questions, but you have to be able to ask them in a concise and lucid way and also need to be able handle the advice being given properly.