Hello all, in order to repreduce an example from a paper, I try to set up a model with the same size of its simulation cell, which is 50.0A × 50.0A × 50.0A. I insert all the molecules by PACKMOL, however this makes too much overlapping and thus leads to a huge total energy in my system even if after minimization. Therefore, I try to use “pair_style soft” command to adjust initial configuration. The total energy did decrease when I used this command but it was still huge when I wanted to switch to actual simulation. I don’t know how to do that. Hope someone can help me. Thanks a lot.
PARM.lammps (1.0 KB)
system.data (1.9 MB)
system.in (1.5 KB)
You have correctly identified the problem: your system has a lot of overlaps and, I may add, a quite high density (but that could be because of the ions).
The approach you chose is reasonable: first relax the structure with a soft potential, then switch to a more accurate one.
To understand what is going on, it is always useful to save a trajectory and to inspect your sample, as it gives a lot of clues. For instance, the potential is too soft and after relaxing the system, there are still a lot of entangled atoms. This is what happens during the minimisation, which is done using the full potential:
A few other comments on your script:
- By using a timestep of 0.1 fs, the system is not changing much. I switched to 1 fs from the beginning.
- You consistently use
reset_timestep 0in all your input scripts. I find it very frustrating, as it makes parsing the output unnecessarily complicated. - A good practice is to print all the contributions to the potential energy, which gives clues on how the sample responds to the various phases of the simulation.
- The strategy I have used is to first lower the density of your sample to give more space to the molecules. Then minimise with the full potential and do a short NVT simulation to dissipate the excess energy, before relaxing the cell vectors with an NPT simulation. The equilibrium density after the above is 1.12 g/cm3.
Here is the revised input script:
system01.in
# LAMMPS input file
variable temp index 298
variable run string system01
variable lbox index 56
# ------------- Initializatiion
units real
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style umbrella
special_bonds dreiding
pair_style soft 1.0
# ------------- read data
read_data system.data
include PARM.lammps
# ------------- set output
dump TRJ all custom 1000 ${run}.dump id xu yu zu
# ------------- soft
neighbor 2 bin
neigh_modify every 1 delay 0 check yes
velocity all create ${temp} 12345
timestep 1
thermo 1000
thermo_style custom step etotal evdwl ecoul elong ebond eangle edihed eimp pe ke temp press vol density
fix 1 all nvt temp ${temp} ${temp} $(100*dt) drag 2
run 30000
# ------------- Lower the density
variable scalex equal $(v_lbox/lx)
fix 2 all deform 10 x scale ${scalex} y scale ${scalex} z scale ${scalex}
run 30000
unfix 1
unfix 2
# ------------- Full potentential + minimisation.
# ------------- It disentangles ions forming pairs.
pair_style lj/cut/coul/long 12
kspace_style pppm 1.0e-4
include PARM.lammps
dump TRJ2 all custom 1 ${run}.min id xu yu zu
minimize .01 .001 1000 100000
undump TRJ2
# ------------- quick equilibration to dissipate the excess energy.
velocity all create ${temp} 12345
fix 1 all nvt temp ${temp} ${temp} $(100*dt) drag 2
run 40000
unfix 1
# ------------- NPT simulation
fix 1 all npt temp ${temp} ${temp} $(100*dt) iso 1 1 $(1000*dt) drag 2
run 50000
unfix 1
write_data ${run}.final
I attach the log and the sample after cooking, so you don’t need to run the simulation again. Courtesy of MaterialX mighty workstation!
As a final remark. Your set-up for this simulation was sound. Is there anyone in your institution to discuss your results? This kind of problems is quite common and takes just a bit of physical reasoning around what is going on during the simulation.
All the best!
Otello
system01.log (50.3 KB)
system01.final.gz (1.1 MB)
1 Like
Thanks for your advice. As a beginner of LAMMPS, there is still something I don’t quite understand. I find you use fix deform command to lower the density, but I think this will change the size of box to be 56A, right? However, in the paper which I want to repreduce, the size of box is 50A, so I don’t know whether this is feasible. Or maybe I can use fix deform command again before formal simulation and after soft process to shrink the box to its original size?
I used fix deform because you cannot expect to get the expected density with a barostat and a soft potential. The density obtained with the (short) NPT equilibration depends on the initial configuration and the full potential. If it’s not what the paper reported, it could be due to a different morphology (e.g. crystalline vs amorphous), to a different pressure, or to a different force field being used in their simulations. Also, never exclude mistakes on peer-reviewed articles: reproducibility is a big issue, especially when raw data, and a clear description of simulation parameters and protocols, are not shared or adequately justified.
Said so, you can try to thermally anneal your sample, that is heat it up to around the melting point and cooling it down at different rates, or use again the deform command to squeeze it to the target density, and then relax it, to see if the morphology changes. Also, you could try generate other starting geometries and repeat the simulation with something like
read_data system.data
pair_style soft 1.0
pair_coeff * * 1.0
variable prefactor equal ramp(1,30)
fix 1 all adapt 1 pair soft a * * v_prefactor
instead of the include you have used. In this way, we are using the same potential for all atoms, and we are increasing the hardness during the first annealing in the NVT step, which should result in a configuration without overlapping atoms.
Again, this is something you should discuss with someone who has some expertise on MD, as the results of computer simulations depend pretty much on those choices. And you are the sole person who is responsible for them.
Now I understand. Thanks for your reply.