Simulating a Solid-Liquid interface with read_data

Hi everyone,

I’m currently trying to simulate a solid-liquid interface to determine the melting temperature of NaCl.

My plan is to have a npt fix act on half of the system to keep it solid at T = 10 K, p = 1 bar and heat the other half to above the critical temperature so that it melts, then run the full system in a NPT ensemble at either slightly above or slightly below the melting temperature and record whether the system is a solid or a liquid and through this way determine the melting temperature.

I have created two data files with the positions and velocities of the two sub-systems. I ran both (the heating/equilibration) separately because the combined system seems to compute much slower.

Solid input file:

/// DEFINITIONS ///

dimension 3
units real
boundary p p p
atom_style charge

/// ATOM GEOMETRY ///

Na atoms

lattice fcc 5.64 origin 0 0 0
region box block 0.0 5.0 0.0 5.0 0.0 5.0
create_box 2 box
create_atoms 1 box

Cl atoms

lattice fcc 5.64 origin 0.5 0.5 0.5
create_atoms 2 box

/// OTHER SETTINGS

mass 1 28.990 # Na mass
mass 2 35.453 # Cl mass
set type 1 charge +1.0 # Na-ion charges
set type 2 charge -1.0 # Cl-ion charges

neigh_modify delay 0 every 1 check yes

/// FORCE FIELDS ///

pair_style born/coul/long 8.0
pair_coeff 1 1 9768.28 0.317 0 24.187 -11.5 # Na-Na interactions
pair_coeff 1 2 28937.74 0.317 0 161.21 -200.0 # Na-Cl interactions
pair_coeff 2 2 80367.73 0.317 0 1669.6 -3353.63 # Cl-Cl interactions

kspace_style ewald 1.0e-5

/// VARIABLES ///

variable dt equal 1.0
variable Tstart equal 10.0

timestep ${dt}

/// ENERGY MINIMIZATION ///

minimize 1.0e-10 1.0e-10 100000 100000

reset_timestep 0

/// SIMULATION SETUP ///

velocity all create {Tstart} 12345 mom yes dist gaussian fix equilibration all npt temp {Tstart} {Tstart} (100.0dt) y 1.0 1.0 $(1000.0dt)

run 10000
write_data Solid.data

Liquid input file:

region box block 0.0 5.0 5.0 10.0 0.0 5.0

variable Tstart equal 10.0
variable Tstop equal 1600.0

/// SIMULATION SETUP ///

velocity all create {Tstart} 12345 mom yes dist gaussian fix melting all npt temp {Tstart} {Tstop} (100.0dt) y 1.0 1.0 $(1000.0dt)

run 150000

unfix melting
fix equilibration all npt temp {Tstop} {Tstop} (100.0*dt) y 1.0 1.0 (1000.0*dt)

run 10000

write_data Liquid.data

As can be seen from the npt-fixes in both files I let the box expand in the y-direction so that the two sides of the interface will have the same area.

I then read both data-files in a third input-file:

/// ATOM GEOMETRY ///

region box block 0.0 5.0 0 14.0 0.0 5.0
read_data test.data group solid
read_data Liquid.data add append shift 0.0 10.5 0.0 group liquid

I’m not sure how much space should be between the two sub systems but probably more than 1.5 Å because below this separation distance the force changes rapidly and causes high accelerations.

I wanted to heat the group solid to T_melt - 50 K in a NPT ensemble and cool group liquid to the same temperature but having two fixes acting on the same simulation box is No Bueno causing the following error:

“ERROR: Must not have multiple fixes change box parameter”

I presume this is because each npt fix will deform the simulation box in different ways causing a conflict. I instead tried to simply program the following:

fix equilibration all npt temp 1024 1024 (100.0*dt) y 1.0 1.0 (1000.0*dt)​ (the anisotropic expansion was an oversight left behind from the previous fixes, not sure if this matters).

Nothing crazy happens when I do this, but how (ill) advised is it to have a different commanded thermostat temperature than the current temperature of the system? Is there a better way to solve this problem?

I also noticed (as can be seen from the line “read_data Liquid.data add append shift 0.0 10.5 0.0 group liquid”) that you need a higher shift value than the space you want because the origin of the second data file is on the opposite side of the interface, I guess.

Best regards,

William Green

trying to glue together two individual subsystem is not going to save you time, but will likely cost you time.

what you are disregarding is that now you have two interfaces in each of your systems that are completely not equilibrated for the conditions of the combined system and that grafting these on top of each other will a) likely cause some shockwaves because of either adding a vacuum (and the system “imploding”) or having suddenly high potential energy and thus the equivalent of “given a kick” to all atoms at the interface, so you need to add dissipative thermostatting to get rid of that, since otherwise it can continue to propagate through your system for extended periods of time and falsify your results and b) now need to wait until the interfaces have “healed” which will happen automatically and with less disruption if you just set up the whole system to below the melting temperature and then just melt half of it with a sufficient temperature delta for the designated liquid part. you can benefit from the extra work that you have done, since you can now estimate the required total volume from your current runs and then don’t even need to use fix npt.

axel

Hi Axel,

Would melting half the system and keeping the other half at constant temperature be possible when having loaded the atoms with read_data? I.e. how do you get around the error

“ERROR: Must not have multiple fixes change box parameter”?

Is it advisable to instead simply assign half the system to one group and the other half to another and skip the read_data stuff altogether?

If group assigning is the preferred method, how can I assign all the atoms in a subspace of the region (say half the y-length of the box)? Is there an easy way to do this rather than having to program it outside lammps?

Best,

William

a) you must not use fix npt multiple times. or any combination of fixes that change the box volume. there is only one box and pressure is a global property and you cannot apply box scaling to just part of the box (you can apply the dilation to part of the system, though, but that also has issues; please see the manual on how to dilate only a part of the system).

b) you can assign groups by region you can also thermalize by region using biasing though compute temp/region. details for both are in the manual.
c) there is no need to use read data for such a simple system. however, you may want to save restart (or data) files at each stage of the system preparation process, so that you don’t have to redo the entire simulation in case you need to make adjustments to a later step.
d) i think you are overly sensitive to the need to use fix npt to minimize error. the error in how well you can determine the point where neither the solid nor the liquid phase is growing or shrinking is so large, that it will dwarf the impact over using a variable cell vs. a constant volume simulation (and that is not even considering the limitations and intrinsic errors of classical potentials and their parameterizations), specifically considering that for your system with half solid, half liquid this is not a well defined process in the first place. since you do have a good estimate for the required volume, it is making your life massively easier to do the process with constant volume.
e) i vaguely recall having given a more detailed outline of how to do a coexistence simulation in the past on the mailing list. you may want to search the mailing list archives.

axel.

Hi Axel,

I do recall you responding to questions about coexistence simulations before, it’s where my question about assigning groups came from. Unfortunately I could not find it when searching the mailing list again.

  1. I did some NPT simulations of either high-temperature solid phases of NaCl or liquid NaCl to get a feel for what the volume should be. I’m somewhat uncomfortable/unsure about doing it this way however, it feels like I’m deciding the density post hoc. If I create the system with a volume that gives the density of a solid, am I not forcing the system to be solid (pressure will increase)?

  2. I tried doing a simulation with NVT ensembles where the volume was the same as a box with half liquid and half solid but after equilibrating (all NVT fixes at each stage) at T_melt - 50 K the system stabilized (at least for the timelength of 100’000 steps after reaching 1024 K for both groups) as a liquid.

I then chose the volume corresponding to the whole simulation box being filled with a solid and was able to achieve a solid phase at the end of the simulation, but the pressure readout was around 6kbar.

Should you increase the volume in a sort of iterative process until you get a reasonable pressure scalar value? Reading melting curves for NaCl it seems like you need several GPa in pressure to significantly change the melting temperature, perhaps this was what you meant in point d).

Here is the most relevant parts of my input file:

variable dt equal 1.0
variable Tstart equal 100.0
variable Tstop equal 800.0
variable Tfinal equal 1500.0
variable Tmelt equal 1024.0
timestep ${dt}

/// SIMULATION SETUP ///

velocity all create {Tstart} 12345 mom yes dist gaussian fix heating all nvt temp {Tstart} {Tstop} (100.0*dt)

thermo 100
thermo_style custom step pxx pyy pzz lx ly lz temp press etotal

dump 1 all custom 10 Interface2.out id xs ys zs
dump 2 solid custom 10 Sol.out id xs ys zs
dump 3 liquid custom 10 Liq.out id xs ys zs
run 100000

unfix heating
fix equilibration all nvt temp {Tstop} {Tstop} $(100.0*dt)
run 10000

unfix equilibration
fix equilibration2 solid nvt temp {Tstop} {Tstop} (100.0*dt) fix melting liquid nvt temp {Tstop} {Tfinal} (100.0*dt)
run 100000

unfix equilibration2
unfix melting

fix heating2 solid nvt temp {Tstop} {Tmelt} (100.0*dt) fix cooling liquid nvt temp {Tfinal} {Tmelt} (100.0*dt)
run 100000

unfix heating2
unfix cooling

fix equilibration3 all nvt temp {Tmelt} {Tmelt} $(100.0*dt)
run 100000

I did still load data-files for the solid and liquid parts so that I could assign the groups that way, but I only ran simulations on the full system. I wanted to run the simulations ASAP, that’s why I chose this perhaps inferior method, but I will read up on the group function in the manual.

3. My last question is about how the atoms moved in the two groups during the simulation. Due to the periodic boundary conditions, some of the atoms from the liquid group “teleported” to the left side (the side of the solid group) and this in turn squeezed the solid group to the middle between two phases of liquid. The system (in the y-axis) looked like this: Liquid - solid - liquid.
My only intuition about this is that it could perhaps speed up the equilibration when both the phases are fixed to the same temperature, but is there a way this could (significantly) negatively affect the accuracy of the results?

Thanks for all the help!

Best,
William

please note that most of the questions you have now are not really about LAMMPS but about how to do your research. that is not something I have time (or interest) to discuss in detail.

Hi Axel,

I do recall you responding to questions about coexistence simulations before, it’s where my question about assigning groups came from. Unfortunately I could not find it when searching the mailing list again.

that means you probably need to improve your search skills. there are likely plenty of matches, since this is a recurring topic, so you may need to look through more than just the first page of results.

  1. I did some NPT simulations of either high-temperature solid phases of NaCl or liquid NaCl to get a feel for what the volume should be. I’m somewhat uncomfortable/unsure about doing it this way however, it feels like I’m deciding the density post hoc. If I create the system with a volume that gives the density of a solid, am I not forcing the system to be solid (pressure will increase)?

since you have already equilibrated a solid system and a liquid system, you can just take the average density and use that. i don’t see why that is a problem?
you have many more ad hoc choices in your system that you are not worrying about, some of which are much more likely to have a significant impact on the accuracy of the results. please also keep in mind, that you are studying a system that is not very compressible, which means that you will see large changes in pressure even from small changes in volume. if you feel this has a significant impact on the results, you need to make tests proving it. please also keep in mind that you have to worry about finite size effects, too; in other words, it may be required to choose a sufficiently large volume at which point the impact of limiting fluctuations due to a fixed volume will be lower.

please note that using NPT isn’t so great either. if you dilate all atoms you are moving the solid atoms beyond where they “want” to be, but if you dilate the liquid atoms only, you have a too large relative movement at the interface between the two groups, i.e. you either create too much of a vacuum or push the liquid atoms too far into the solid atoms. if you don’t dilate the atoms at all but only adjust the box, you have the same thing but only at double the magnitude at the box boundary. using a fixed volume will eliminate this undesired side effect of NPT (or NPH).

  1. I tried doing a simulation with NVT ensembles where the volume was the same as a box with half liquid and half solid but after equilibrating (all NVT fixes at each stage) at T_melt - 50 K the system stabilized (at least for the timelength of 100’000 steps after reaching 1024 K for both groups) as a liquid.

then you were using the wrong temperature. the whole point of a coexistence simulation is to find the temperature where the size of these parts of the system do not change (significantly).

I then chose the volume corresponding to the whole simulation box being filled with a solid and was able to achieve a solid phase at the end of the simulation, but the pressure readout was around 6kbar.

Should you increase the volume in a sort of iterative process until you get a reasonable pressure scalar value? Reading melting curves for NaCl it seems like you need several GPa in pressure to significantly change the melting temperature, perhaps this was what you meant in point d).

Here is the most relevant parts of my input file:

as mentioned about. i don’t have the time to review the science, nor is this a topic i have recently worked on, so to give any meaningful advice specific to this kind of system, i would have to spend a significant amount of time to read up on it and make tests to reconfirm my previous experience. this is obviously beyond the scope of this mailing list and something that i don’t have time to do. so you have work out how to do proper science with your adviser.

[…]

3. My last question is about how the atoms moved in the two groups during the simulation. Due to the periodic boundary conditions, some of the atoms from the liquid group “teleported” to the left side (the side of the solid group) and this in turn squeezed the solid group to the middle between two phases of liquid. The system (in the y-axis) looked like this: Liquid - solid - liquid.
My only intuition about this is that it could perhaps speed up the equilibration when both the phases are fixed to the same temperature, but is there a way this could (significantly) negatively affect the accuracy of the results?

nope. a 3d periodic system is theoretically invariant toward translation. normally, that would require removing the center of mass momentum after equilibration. using fix npt will increase the likelihood of forming a center of mass drift due to floating point rounding inconsistencies, if the origin of the cell is badly chosen. with a part solid, part liquid system, this gets worse. for as long as this is small, it is nothing to worry about, and it can be easily “fixed” after the fact by retranslating the system so that the center of mass of both sections remain in place. please note that with a periodic system you will always have a liquid-solid-liquid-solid-liquid-solid-liquid-solide-… system and the choice of origin is arbitrary. if you feel it is too much of a burden, you can add a soft tether with fix spring to the solid phase to keep it from drifting, or use fix momentum to remove any drift or use fix recenter to hide the drift (but beware of the “flying icecube syndrome”, which may be hidden by the latter command. in short this is a none-issue and has no impact on the performance of the calculation.

that is about everything i can say about this topic. don’t expect me to answer to anything else outside of technical questions that are concerning LAMMPS directly.

axel.

​Hi Axel,

Thank you for answering the question I had that concerned Lammps. Your comments about the more general issues were helpful too, even though as you said those questions are outside the scope of the mailing list.

Best,

William Green

Dear Lammps users,
Recently I tried to heat up some atoms in a circle region with a temperature gradient (there is a highest value at the center of the circle) from room temperature and then cool down the whole region to a room temperature.
I tried:

compute new powder temp/region first
compute ke1 all ke/atom
compute pe1 all pe/atom
variable temp1 atom c_ke1/(1.5*8.617343e-5)
compute msd1 first msd
compute stress1 all stress/atom NULL

please note that temperature is a bulk property and only well defined in the limit of infinitely large systems.
when computing temperature from kinetic energy, we assume equipartitioning and disregard (normal) fluctuations.
if you look at the most extreme example, a single harmonic oscillator, you have the “temperature” switch between 2x the average and 0.
so expecting well defined and smooth temperature gradients for local structures is neglecting the underlying physics and specifically statistical mechanics and thermodynamics.

that said, if your goal is to thermalize to a given target temperature in a specific region, you can use a temperature bias using compute temp/region, which will aim to set the temperature in the given region and not deposit a predefined amount of kinetic energy. details for how to use temperature biases are in the manual.

axel.

Thanks very much for your excellent advice.

Also note that if atoms/particles are moving in and out of your region,
e.g. for a liquid, then you need to use a dynamic group to heat
(or thermostat) only particles in the region. And likewise to
measure the temp of only atoms in the region.

Steve