Glass-forming liquid

LAMMPS LAMMPS Beginners
1

Hi everyone,
I am going to simulate a glass-forming liquid using a model described in the article:
https://www.science.org/doi/full/10.1126/sciadv.aba8766
(It models a glass-forming liquid using a mixture of two types of particles (80:20) interacting with each other through a smoothed LJ potential). This model shows the transition around 0.35 temperature in the reduced units. However, my simulation finds this material in the glass phase even at 0.466 temperature!
I have been struggling with this issue for a couple of weeks! and tried to modify every parameter that I thought could play a role in this behavior. I examined other thermostats, Langevin, Berendsen, …, I also test rescale velocity method and fix nvk which were mentioned in references of the article or the related ones. All methods can keep the system at the desired temperature, however, the system still goes into the glass phase at higher temperatures than the critical temperature.

I provide a simple version of my Lammps code here. (this version applies lj/cut as pair_style potential, but my original code calls a table letting the particles interact with each other through the smoothed lj/cut. Also, I eliminated the preparation of the initial configuration section).
I would be really thankful if you could guide me with this issue.

variable rho equal 1.2
variable Natom equal 1000
variable L equal ({Natom}/{rho})^(1/3)

units lj
dimension 3
boundary p p p
atom_style atomic

region SimulationBox block 0 {L} 0 {L} 0 ${L}
create_box 2 SimulationBox
create_atoms 1 random 800 737392 NULL overlap 0.5 maxtry 10000
create_atoms 2 random 200 647483 NULL overlap 0.5 maxtry 10000

neighbor 0.3 multi
neigh_modify every 1 delay 0

mass * 1.0

velocity all create 0.466 837376 dist gaussian mom yes

thermo 5000
timestep 0.01

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2

minimize 10e-7 10e-7 100000 100000

fix fix1 all nve
run 50000
unfix fix1

fix fix1 all nvt temp 0.466 0.466 1.0
run 100000
unfix fix1

fix fix1 all nvt temp 0.466 0.466 1.0
dump dump1 all xyz 5 dump.xyz
run 100000
undump dump1
unfix fix1

2

Sorry, but this is not a question about LAMMPS but rather about how to do your research project. It is thus off-topic for this category and rather a topic to discuss with your adviser and your tutors and collaborators . Please see the forum guidelines for the LAMMPS categories: Please Read This First: Guidelines and Suggestions for posting LAMMPS questions
Please also consider the notes about formatting your messages in future posts.

3

Dear Axel,
sorry for the inconvenience. Actually, I was not sure if I could post my question here. however, I did it because I thought it might be a very familiar issue for experienced people in Lammps.
(Indeed, I have been searching for help for a while but I found none).

4

If you carefully observe the ebb and flow of the forum, you should have noticed that most topics are created by people with rather limited experience that mostly have “technical problems”, either very fundamental or within a given functionality of LAMMPS. You should also notice that the (rather few) people that respond to posts, can for the most part only help with these technical problems, since they are unfamiliar with the research topic. Considering the wide range of functionality in LAMMPS chances are extremely slim that you come across someone that has done research on the same topic.

As a general rule, what you do not find here much are experienced practitioners in simulations. They usually do not need to get help in running LAMMPS and they are usually busy hustling for research grants and advising their own students and postdocs, that they are not out on the web advising other people’s students or postdocs. Since research is in may ways a business these days, that would be similar to a PepsiCo executive helping people working for CocaCola. This is a rather unfortunate affair. So your only chance is to get help from people in your “company”, and frankly, that is their job.

5

I didnt really read that paper you linked carefully nor I am familiar with some terminology the authors are using, but by judging by the part shown below, I would say the potential form you are using is not the same as the one they used (which may be the one reason why your results are different).

Maybe you can search in the manual ways to implement a “customized” potential form of your choice (I dont know if the potential shown below exists in LAMMPS).

image
image901Ă—439 31.4 KB

1 Like
6

Thank you Axel for your time and attention to my problem.

7

Thank you ceciliaalvares for your time to help. However, the potential you have mentioned is related to another part of the simulation where a unique particle is added to the simulation.

8

ahhh, makes sense, my bad :"D

Well, assuming that the potential form is correct for what you are trying to do then, what I can see from your script that can be a bit weird is the minimization you are doing after having assigned velocities to the particles in your system (e.g. of this problem: Energy minimisation (CG) at various temperatures). Maybe you can also check the temperature evolution during the nve portion of your system (PS: I dont see why you would need to do some dynamics in the nve ensemble nor an energy minimization for glass transition, but indeed I am not familiar with the topic, so as Axel said this is probably seen better with your supervisor or the team you work with). I think for glass transition you start at a high temperature (at which the liquid phase is the stable one) and then you drop it very quickly or do a ramp where the temperature changes very fast, no? So maybe you can equilibrate your system in the liquid phase first (again, to be seen with your supervisor).

You can also save data files (see the write_data command page) and dump file before the last part of the simulation to try visualize what’s happening and carry out a deeper investigation yourself.

If you dont have anyone around to assist you with the project, you can try asking questions on the Science Talk category about “how people working with glass transition simulate it” or try to see if someone working with this happens to make input scripts available on github or something for you to try to learn from it.

9

I think this one is it: pair_style lj/expand command — LAMMPS documentation
The \frac{1}{4} is just a shift same as pair_modify shift yes.

1 Like
10

Hi @Fatima,

@ceciliaalvares was wrong about the potential but right in the fact that you are not using the smooth potential of the paper.

The formula they use is
U(r) = 4\epsilon[(\frac{\sigma}{r})^{12} -(\frac{\sigma}{r})^6] + 4 \epsilon[\sigma_0+c_2(\frac{r}{\sigma})^2 ], see equation 6.

With each terms depending on the species interacting. This differs significantly from the lj/cut formula.

A don’t know if this pair interaction is implemented in current LAMMPS stable or dev release. But as suggested using pair style table should be rather straightforward.

11

Dear ceciliaalvares, Thank you for your comments. Minimization changes the coordinates of the particles to reach a configuration with the lowest energy. So velocities have no effect on it. no? I use fix nve to remove the initially assigned velocities. It is a step that is written in the article. If fix nve is used alongside the Langevin thermostat, it can set the system on the desired temperature too.
About glassy systems, you are right. They start from a high temperature and quench it. My problem is that my system is set at a temperature higher than glass transition, but it behaves like glass!
thank you for your suggestion. I will ask my question.

12

Thank you Germain for your reply.
As I mentioned in the first post, I have provided a table for the smoothed LJ interactions taking into account the extra two terms satisfying the continuity of the potential and force. The code here is a simple version of it. My problem is still there for both original and smoothed LJ potential.

13

Sorry for reading your initial posting too quickly. I understood that you were trying to find out your problem using this script.

The problem is then that it makes very hard for us to help you since it is difficult to know what you exactly do with regard to the paper. Looking more carefully into your script, I see major differences with the protocol of the paper. For instance it is said that the NVE relaxation was done using the soft pair_style and it is not said that any initial velocities were provided to the atoms. This is absent of your script while stated to be a crucial step of the paper.

You could switch with the following commands*:

pair_style soft 2.5
pair_coeff * * 1

fix fix1 all nve
run 50000
unfix fix1

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2

However it is impossible to determine if you did this or not since the minimal working example provided apparently differs too significantly compared to what you actually do. Sorry for not providing better insight.

On your behalf, I have to say that the protocol provided by the paper is very… short… to say the least. For example the A constant value of the soft potential is missing.

It is hard to help you further here and as stated by cecilia, you might get better insight from your local co-workers or reading more about how simulation of glassy materials is done in general.

2 Likes
14

I provided a simple version to make it straightforward. I should say that I could not upload the table file too. I do not know if there is another way to share it with you?

And about the major differences you pointed out, you are right. In the first post, I mentioned that the preparation of the initial configuration section is eliminated. Not only the relaxation part using the soft pair_style but also shearing by nvt/sllod. It is because as far as check these steps can not modify the glassy behavior of the system that I observe at high temperatures, 0.466 for example. On the other hand, as you noticed the protocol explained in the paper is very short. A constant, rc constant, and shearing parameters, all of them are unknown.
I should let you know that a while ago I could make contact with one of the authors of the article and asked about the parameter values. He explained that the initial preparation steps are not necessary. Instead after assigning the initial random velocities (according to the article) a long-time run at high temperature (they mentioned 0.466 in their article) can relax the system. Unfortunately, I could not reach him anymore. Generally, no one is around at the moment.

15

Sorry, I am not exactly at my sharpest day and unlike Germain I am being a bit too lazy or incompetent to properly read the stuff from the paper :“”")

I would try doing the preparation stage (with the soft potential they mentioned) to see if it changes something. Maybe the author you contacted was not the same person who wrote the sentence about how crucial the preparation stage is and it is possible that the latter person knows better than the former - who knows… besides, it’s not like you have many other options left.

Minimization changes the coordinates of the particles to reach a configuration with the lowest energy. So velocities have no effect on it. no?

I dont know the specifics of what happens in the energy minimization protocol upon assigning kinetic energy to the atoms, but apparently there is a problem with proper Ep convergence since the information can be found this stated in different posts on the forum. If that’s the case why not try to do it before assigning the velocities? But also there is no energy minimization in the methodology, right? (not that it should make a difference to have it I suppose)

Did you check if the temperature remains at your desired target one during the dynamics in the NVE ensemble? It could be that the phase transformation is happening there. This would support the idea of the “preparation of the system” being important.

16

I use fix nve to remove the initially assigned velocities.

Could be that my brain is still being too slow, but I am not seeing how doing some dynamics in the nve ensemble is going to remove the initially assigned velocities… (and even if it is going to do that, why assign velocities in the first place if the goal is to get rid of it?)

nazare
nazare702Ă—460 50.3 KB

Unless you are having a phase transition that is consuming all the existing kinetic energy of your system…

EDIT: I understand now: you mean that it will reduce, not remove the velocities. This happens as a consequence of you not initializing the system with the due potential energy to keep the Ek when running the dynamics in the nve ensemble. Note that it could also increase the kinetic energy of your system depending on your initial value of Ep, so it doesnt mean it will necessarily decrease the velocities.

If fix nve is used alongside the Langevin thermostat, it can set the system on the desired temperature too.

But there is no fix langevin or anything coupled to the fix nve, so I suppose it is indeed a simulation in the nve ensemble in your case.

17

Have you tried starting your simulation using the setup you showed us with an initial microstate that maintains a correspondence with the occuring ones at a thermodynamic state of (N = your N, V = your volume, T = 0.466)? This means that it would also need to have somewhat the proper corresponding internal energy associated to this thermodynamic state (with possible low fluctuations).

My point is: if your initial microstate is not one coming from that thermodynamic state (and is actually far from it), maybe your temperature is changing during the first stage of your simulation (corresponding to the fix nve), possibly in the direction of diminishment, as Ep and Ek interconvert. And this could actually be the reason why you are observing an amorphization. It is as if the amorphization happens in the correct range of temperature and during the nve part but you are observing it at “higher temperature” because you are saving data only corresponding to the later part of your simulation.

Try removing entirely this part below:

minimize 10e-7 10e-7 100000 100000

fix fix1 all nve
run 50000
unfix fix1

and run the simulation to see if it will be the amorphous or liquid structure. If indeed that “preparation stage” is not necessary, my bet is that it will be the liquid state :slight_smile:

18

Thank you Cecilia for the series of your comments. I am checking them one by one.

About the minimization, It just takes action after the pair_style definition. Before that, nothing happens. It is an essential step because it aparts overlapped particles. Indeed, without this step, the system explodes with or without “soft pair_style” at the beginning of the simulation.

I constantly monitor the temperature during the simulation. when I use nve ensemble before applying nvt ensemble, The temperature is lower than the desired 0.466 whether the interaction is “soft” or “LJ”. It makes sense because of the minimization we did.
Using a thermostat like the Langevin thermostat and nve at the same time can set the temperature to the desired value. This is different from what is found in the protocol. They might change the “soft” pair_style magnitude to achieve the desired temperature or leave it to the next step, nvt ensemble, to get the specific temperature.

19

I am not sure if I get your point. Initial velocities are random and nve dynamics let them adjust their velocities according to the force acting on them. Maybe you mean this step is not essential because the next nvt integration will adjust them. Honestly, I am not sure what is the role of this part.

20

Which observables do you monitor to confirm that the system is in the glass phase?

Additionally, have you tried simulating 4000 particles (as per the paper) instead of 1000 particles? Sometimes finite size effects can make a big difference.