How to simulate diffusion of hydrogen molecules within a simulation box?

Previously, I have successfully simulated carbon dioxide molecules within a simulation box, by assuming carbon dioxide as a rigid model.

Now, I want to simulate hydrogen molecules, but there are errors. The errors include: “Fix temp/berendsen period must be > 0.0” and “Lost atoms: original %ld current %ld”. I have read through some journals, and they stated that hydrogen molecules should have a center bond charge, which means that there would be three charges for a hydrogen molecule (2 for the atoms, 1 for the center). The problem arose when I assumed the center bond charge as an atom in the hydrogen molecule file. I am not sure whether this is the problem contributing to the errors, but I know it is weird to start with this assumption. I have also come across some forums, and some of them discussed that LAMMPS isn’t able to simulate center bond charge. However, if I need to simulate the hydrogen molecules using LAMMPS, could anyone here suggest to me how should I proceed? Does anyone in the community have successfully simulated hydrogen molecules before? The hydrogen lammps input file and its molecule file are as follows:



Thank you.

Please do not copy past images here, but instead the raw inputs files with proper formating.

How did you choose the constants for the bond and angle coefficients? A value of 10000 for both quantities seem completely random. Also note that the minimum timestep you can use in a simulation is set by the vibration frequency of the bonded atoms, so using random large bonding coefficients will likely make your simulation fail.

Thats not the only big problem with your script btw, you also have 2 thermostats acting on the same atoms.

Simon

Dear Simon, I am a new user for LAMMPS, I can’t post any files in the forum just yet, which was why I chose to post the screenshot of the scripts. I am currently using a rigid hydrogen model, that has two LJ-sites models with three partial charges. One of the charges is in the center of the bond between the two hydrogen atoms.

Since the hydrogen model that I set was rigid, the bond coefficient I set was 0. However, there was an error. So, I tried different values, such as 10, 100, and 10000, and all of them showed errors. I even tried to put negative values, but the error persisted. If I put zero or values under zero, the error would be “Lost Atoms”; if I replace the values with values higher than zero, the error would be “ERROR: Computed temperature for fix temp/berendsen cannot be 0.0 (src/fix_temp_berendsen.cpp:142)”.

I sincerely hope that you can assist me with my simulation because I am currently very new to molecular simulation. The main problem that I want to solve is to simulate a hydrogen model that has a charge at the center bond. Second, I hope that you would clarify to me why 2 thermostats could not act on the same atom. At first, I thought of trying to use NPT ensemble, and use the Berendsen to control the temperature. Is it redundant to use Berendsen to control the temperature having the NPT ensemble with a thermostat? Thank you.

Regards.

You can use the preformated text option to write LAMMPS script, e.g.:

units reals

I don’t see any command that would make the molecule rigid in your script.

All the issues that you are having suggest that you need to take a step back and go through tutorials before attempting to make this simulation, otherwise you will keep having new issues again and again even after solving the ones you have right now. You may also want to learn about the fundamentals of MD, for instance by reading “Understanding molecular simulation” by D. Frenkel. Answering these issues is not the goal of this forum I am afraid.

Best,
Simon

This is strange since that kind of model would be effectively the same as a \mathrm{CO}_2 molecule, so you should be able to adapt from your previous simulations by using the same basic approach, but just replace the O atoms with your H atoms and the C atom with the dummy charge point.

LAMMPS doesn’t know (or care) about the chemistry, it just implements the model you define by looking at points, distances, charges, masses, and then computes forces and processed those with the fix commands to define.

That statement is not entirely accurate. This statement also is not very useful since you are not providing any references. Remember, we are trying to do science here and no statement like that can be trusted unless you provide proof or a reference. How can we know that you were not misreading what you saw?

You were right. I will go through the lammps tutorial first before coding my script. Thanks for your advice.

Regards.

I will review back my script. Thanks for your advice.

Sorry for the misjudgment. I was then searching for a solution regarding the lost atoms issue of my hydrogen model, and I found someone saying this, which I could have been totally mistaken. Will not judge facts without evidence ever again. Thanks for the advice.

Regards.

Without looking at this thread in too much detail I would recommend looking at fix gcmc to see whether it’s appropriate for whatever problem you’re trying to solve