I am trying to model iron (III) oxide surfaces to be used in a simulation of a lubrication system.
I have tried the following approach: I used the crystallography open database to get the cif file for the crystal and then used “Atomsk” to generate a larger cell and then cut the cell in half horizontally to have both bottom and top layers.
Following that I used the outputs from “Atomsk” to generate the initial coordinates of the surfaces using “Packmol” and that is where the problem became much more evident: naturally the crystal has a rhombohedral shape and due to this shape, there would be areas that are uncovered by it in the simulation box.
I would attach the xyz file to better showcase what I am referring to but I lack the permission to do so as a new user on the forum.
As such, how would I resolve this “issue” using lammps, or should I use any exterior software, if so then what are your recommendations? I hope it is not a problem to be asking about using exterior software as this is a lammps forum, if it is, my apologies.
Isnt this something you can address by using a triclinic simulation domain? It would suffice to prepare a LAMMPS data file with the proper dimensions and angles (see read_data command — LAMMPS documentation for more details on how to set up a triclicnic simulation domain)
I don’t see why you would need to use packmol here. To the best of my knowledge Atomsk can already output atomic structures in LAMMPS data file format.
Please also note that you can create such structures which are based on simple lattices directly in LAMMPS use the lattice custom and create_atoms commands. However, in all cases some knowledge in basic crystallography is needed, so if you don’t have much experience there, you may want to find a colleague that has it and is willing to help you and explain things to you.
Please note that metal oxide surfaces are a tricky business because of the polarity of the atoms/ions. This will make it very important (unlike with mostly neutral metals) how exactly you cut the surface, and you will have to deal with reconstruction and avoiding net charged systems. Also, you will likely need to look for rather complex (and thus slow to simulate) potentials with charge equilibration to model this kind of system.
Thank you very much for your response, let me provide more information about my reasoning for using packmol:
The system that I want to simulate is a lubrication system where lubricant molecules (hexadecane and or fatty acids representative of a vegetable oil) will be sandwiched between 2 iron oxide surfaces in order to simulate shearing of the lubricant layers. In subsequent simulations, the aim is to add different types of nanomaterials to the system to hopefully better visualize and understand their mechanisms of interaction with the surfaces and the lubricant molecules, because the experimental studies that I did yielded interesting results and I hope to understand the mechanisms behind them. As such, by looking through the literature, I was able to find some studies that I want to recreate and adjust to fit the purpose of my work. The studies in question:
Molecular dynamics simulation of surface energy and ZDDP effects on friction in nano-scale lubricated contacts: Redirecting
Structure and Friction of Stearic Acid and Oleic Acid Films Adsorbed
on Iron Oxide Surfaces in Squalane: https://doi.org/10.1021/la404024v
It is for this purpose that I wanted to use packmol to generate the initial configuration of the system both with the lubricant molecules and surfaces. Concerning the lubricant molecules, it has so far been “smooth sailing” but my issue has mostly been with the surfaces themselves.
Thank you for your note, I was planning on using the LJ parameters that were tested by the researchers in the papers that I previously mentioned. Would that be the correct approach here?
In your opinion, would it be feasible and scientifically correct to substitute the iron oxide surfaces with a simpler alpha iron surface?
You can create two different data files for the two parts of the system and have the same box dimensions and just leave space for the other part. Then you can add the second data file after reading the first. Another option would be to put the definition of the molecules into a molecule file and then create both, surface atoms and molecules with create_atoms commands. Both approaches have limitations and require careful attention to details of the force field and atom and topology type requirements. Neither of those, nor what you are doing are suitable for beginners in LAMMPS, but require experience and sufficient understanding of how LAMMPS processes types and topologies.
This is not my area of research and thus I am not qualified to give an authoritative answer, but using only Lennard-Jones without charges would be highly suspicious to create bogus results.
Again, this is not my area of research and I am not your adviser and furthermore this is not a question about LAMMPS but about your research. You need to discuss this with the people that have a vested interest in your research. All I can say is that alpha iron and iron(III)oxide are very different materials.
My apologies for not being clear enough, the studies in question provide partial charges as well for the Fe and O atoms, 0.771 and -0.514 respectively.
Thank you very much for your responses, there is a lot to consider!
From a practical point of view (been there, done that), I think the suggestion of @akohlmey of using custom cells inside LAMMPS is way easier to implement in your workflow than using Packmol or other external packages.
This would allow you to define regions in terms of cell length instead of absolute length and might make things a bit easier to navigate.
Consider the following:
lattice fcc 5.0
variable nx equal 10
variable ny equal 10
region full block 0 ${nx} 0 ${ny} 0 20
region bot block 0 ${nx} 0 ${ny} 0 2
region mid block 0 ${nx} 0 ${ny} 2 18
region top block 0 ${nx} 0 ${ny} 18 20
This allows you to create a box using the full region and only use create_atoms command in bot and top for your surface, and the mid region to insert molecules using molecule templates.[1] Here this would results in a box of units 50\times50\times100\ \mathrm{A}^3 box with the 100 axis along z having FCC crystal only in the bottom 10 and top 10 angstrom.[2]
You can replace the fcc lattice with a custom one as shown in the lattice command documentation, and the basis vector and atomic positions are rather straight forward to derive from a .cif file, see the How to triclinic page. Once set, you can also easily rotate the cell to adjust it to your liking.
As a side note, I would recommend to avoid using several pre-processing software especially if you are a beginner in using a simulation package, since it will make analyzing the origin of errors and/or simulation results more confusing. The use of lattice commands is often overlooked but often more straightforward than a combination of several “jacknife” programs.
By default LAMMPS defines the distances used in several commands through crystalline lattices if the units box option is not used. ↩︎
Also assuming you’re not using lj units in which case lattice commands behaves differently. ↩︎
Than you very much for your suggestions @Germain@akohlmey@ceciliaalvares, they were extremely helpful. I prepared a script to model a single unit cell first to verify if my workflow is correct before upscaling it to the larger system that I described before. This is the script:
However when I ran it, I got the following error ERROR: Unknown lattice keyword: & (src/lattice.cpp:224), even though I think I wrote it like how it appeared in the documentation, and in a previous script that I tested before ran fine when using the exact same “&”. What am I doing wrong? Also, is my usage of the region and create_box commands correct here?
Yes, thank you for the advice! it’s just a residual habit from before as in the beginning I started using moltemplate and Packmol to prepare much simpler systems than this, but when I made this switch to working directly with lammps I should have made the appropriate adjustments.