Mass density profile

Hello friends,
I used GCMC to investigate gas adsorption on my adsorbate. I have the final structure, including adsorbate and gas molecules. Now, I want to plot a gas density profile on the length of the adsorbate, and I want to know how I do it.
I searched in the tutorial and found Fix ave/chunk, but this is not what I want! I think this command should be used along with the GCMC investigation, and the result of this command gives us density (1/volume) vs. distance, but I want density with g/cm3 unit or mass density.
Do I need to use other software like VMD or OVITO for this purpose? Or should I do the adsorption process with GCMC + another command?
I would be grateful if someone could help me.
Best regards,

I think this command should be used along with the GCMC investigation

Why exactly do you think that? It’s not mentioned in the documentation.

As for the mass density, did you see the density/mass keyword in fix ave/chunk?

Please re-read the documentation more carefully.

Dear Michael,

Thank you for your answer. You’re right. I should read the documentation more carefully.
But, I hope there is another way instead of doing simulations again!
Doing the simulations again is time-consuming since I have already done them!

Assuming you have a dump file, you can use the read_dump command to compute the density profile without running the simulation completely. You could also use the dump file in visualization software like VMD or Ovito, though I don’t know the exact steps for that.

1 Like

Two more comments:

  • LAMMPS has the rerun command to do additional analysis on an existing trajectory
  • it is best practice to do test runs of projects with very small systems and shortened trajectories to check if an input or model or method works well enough and provides all data needed. You just learned the hard way why this is so
1 Like

Thank you, Axel and Michael, for answering my question.

As you said, I reran and used Fix ave/chunk. But, I have a problem: my simulation box is not rectangular; there are more or less than 90 degrees angles. So, when I use the Fix ave/chunk, I receive an error saying use reduced; when I use reduced, I get results between a range of 0-1! How can I report the results? Versus the range of 0-1 or the number of chunks?
On the other hand, I used the box instead of the reduced and deleted the row of angles in the data file. I get results at the length of my simulation box, but I don’t know if this solution is correct because the code didn’t run with the angles row in the data file!
I would appreciate it if you could give me advice.
Best regards,

You need to be more specific.

What was the error exactly? What exact commands did you add/change to fix it?

What do you mean you “deleted the rows of angles”? Do you mean the box tilt factors? Or do you mean 1-3 bonded interactions defined by the angle_style and angle_coeff commands? You should use the terminology defined in the documentation.

Dear Michael,

when I used this command:
compute densityco2 co2 chunk/atom bin/1d x -2.447001 38.065

I got this error: Compute chunk/atom for triclinic boxes requires units reduced

My point about “deleting the row of angles” is that I deleted the fourth row at the following in my data file:

-2.447001 38.065000 xlo xhi
5.671652 46.550348 ylo yhi
9.269484 54.355516 zlo zhi
-9.957595 12.982013 9.429407 xy xz yz <----- this row

and used this command: compute densityco2 co2 chunk/atom bin/1d x lower 0.5 units box

You dont want to do that. Atom coordinates have to be properly converted from triclinic to orthorhombic.

To properly convert triclinic coordinates into orthorhombic, there are some external tools (I did not tests them, but I suppose they work), such as : Trajectory transformations — lipyphilic.transformations — lipyphilic 0.10.0 documentation

1 Like

Dear Simon Gravelle,

Thank you for your suggestion. I’ll test them.

That’s not really going to work, as your fix is being rerun over the dump file(s), which includes the box information.

The best way forward is to use the reduced units. Take note of the meaning of the delta keyword, etc. with reduced units. Also note the shape of the bins for triclinic boxes. The 0 bin is the lowest bin in whichever dimension, and the 1 is the highest. If you want, you can scale the bin edges by the box width in post-processing, but again, use common sense. The interpretation of the data is up to you.

1 Like

I agree, this is not gonna work if the modification is applied only on the .data file and then the rerun command is used. But with external tool like MDAnalysis, it is possible to import a dump file, apply a transformation (like convert triclinic coordinates into orthorhombic), and then print a new dump file. In a way, this is the same as using the ‘trjconv’ command of gromacs. I’m not really familiar with the rerun command, but I don’t see why it wouldn’t work.

1 Like