I am trying the Ar G-K example at 8.3.6. Calculate viscosity — LAMMPS documentation and everything runs fine, but the system is a solid at these conditions in my run (obvious from zero diffusion and looking at trajectory video in vmd). Everything though seems to run fine. Is there a problem with this example, or what? Thanks for any help!
Sincerely,
Ron Cohen
Since nobody has responded to this so far, I suggest you post your observation as a bug report issue in the LAMMPS GitHub repo. That will give it more visibility for the people that wrote or edited that segment of the documentation and (most importantly), it will not vanish unless it is explicitly dismissed.
Good point. That example uses the nominal melting temperature of argon. I found that I had to raise the initial temperature to 250K to get the structure to melt. Even then, the system did not remain in the liquid state below 200K. I could update the example, but it might be better to redirect readers to the official set of examples in ./examples/VISCOSITY
Thank you so much for the response! Still there is a question here. Why isn’t it molten? I looked up the literature for L-J for Argon and found in THE JOURNAL OF CHEMICAL PHYSICS 143, 234504 (2015)The Lennard-Jones melting line and isomorphism D. M. Heyes and A. C. Brańka that exactly these parameters should be in the liquid at the triple point. So I think the problem is that it is started as a crystal and even though it thermodynamically should be a liquid it does not melt spontaneously . It worries me though that one has to superheat so much ( I found the same thing). So I wanted to try starting with random atomic positions, but get an error when I try:
dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 6 0 6 0 6
create_box 1 box
#create_atoms 1 box
create_atoms 1 random 864 23876921 NULL overlap 1.5
I get: LAMMPS (17 Feb 2022)
using 4 OpenMP thread(s) per MPI task
Lattice spacing in x,y,z = 5.376 5.376 5.376
Created orthogonal box = (0 0 0) to (32.256 32.256 32.256)
1 by 1 by 1 MPI processor grid
ERROR: Illegal create_atoms command (src/create_atoms.cpp:224)
Last command: create_atoms 1 random 864 23876921 NULL overlap 1.5
which I am sure is my ignorance on how to set this up. Any help on this would be greatly appreciated.
So to summarize, I would like to be able to run Ar liquid at the triple point, just like the example tries to do, but successfully!
Sincerely,
Ron Cohen
Yes, the starting solid phase remains metastable well above the true melting point. There are lots of reasons for this (finite size, finite time sampling, etc.). It’s not a good idea to use random positions for a dense fluid. Better to start from a solid and melt it by heating to a sufficiently high temperature, then slowly cooling.
Thank you! Yes I superheated to 300K and then immediately to the conditions in the example, but it crystallized and is not a liquid. The diffusivity is zero. I used:
#melt first
velocity all create 300 102486 mom yes rot yes dist gaussian
#fix NVT all nvt temp $T $T 10 drag 0.2
fix NVT all nvt temp 300 300 10 drag 0.2
run 8000
unfix NVT
reset_timestep 0
fix NVT all nvt temp $T $T 10 drag 0.2
run 8000
So it seems it is crystalline at these conditions in LAMMPS, though the literature says liquid. I will run a larger systen too. I would appreciate any suggestions.
Sincerely,
Ron
I tried now 10x10x10 and melted it at 300K then cooled instantly to 86.4956 K and it is at best amorphous. There is no diffusion even in a million step run for 10x10x10. So its seems LAMMPS L-J doesn’t agree with literature on this? The cutoff here is 13 which is similar if not identical. Could long range corrects make such a difference? It is very puzzling to me,
Sincerely,
Ron
It’s not puzzling. There are many technical details that must be considered when reproducing specific literature results. For example density. Also, detailed settings for the pair_style are very important, tail corrections, etc. Assuming the literature source is correct, there is no reason to expect that LAMMPS is incapable of reproducing the result. It’s just a question of taking the time to eliminate the various differences in simulation protocol. LAMMPS has many degrees of freedom in addition to the 6N positions and velocities.
Thank you. I ran for longer at 300K and at 86K (10000 steps eachs) and now it seems to be a liquid. Here is the r^2 graph versus time.
It gives D=8.15e-10 m^2/s if I have my units correct. I not change anything else. So the top of the fixed example file now shows:
# LAMMPS input script for viscosity of liquid Ar
units real
variable T equal 86.4956
variable V equal vol
variable dt equal 4.0
variable p equal 1000 # correlation length
variable s equal 5 # sample interval
variable d equal $p*$s # dump interval
# convert from LAMMPS real units to SI
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m}
# setup problem
dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.2381 3.405
timestep ${dt}
thermo $d
# equilibration and thermalization
#melt first
velocity all create 300 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp 300 300 10 drag 0.2
run 10000
unfix NVT
reset_timestep 0
fix NVT all nvt temp $T $T 10 drag 0.2
dump myDump1 all custom 10 dumpT${T}.NVT.lammpstrj id type xu yu zu vx vy vz
dump_modify myDump1 sort id
run 10000
I understand caveat emptor and all that, but it still seems it would be good for examples on the official LAMMPS web site to work as given. It would be good for there to be the expected answer there too. Thank you again.
Sincerely,
Ron
Good! Thanks for pointing out the original problem of non-liquid argon, which was not acceptable. We have updated the doc page with a script that melts the argon at 250 K and samples viscosity at 200 K. That is good enough for a demonstration. Collected small fixes by akohlmey · Pull Request #3327 · lammps/lammps · GitHub