Thermodynamic output of simulation on two independant regions

LAMMPS LAMMPS Mailing List Mirror
1

Hi! I’m trying to simulate melting through the interface method, where half of my simulation is liquid and the other half is solid, them being put next to one another should spread the liquid state through the solid half. However I’m unsure how to understand LAMMPS’s output or what to expect from it, and was hoping you could be of assistance! :slight_smile:

Below is my LAMMPS input script, my logic was to create a huge box and split it in two, I’d bring both to 1400K using NPT temperature ramps, going 50K to 1400K for the solid box and 4200K to 1400K for the liquid box, then keeping them at 1400K.

My fix 5 has the purpose of printing the thermodynamic results found in the log file, however it seems to only be printing the thermodynamic results of the “solid” box, as it only contains temperatures between 50K and 1400K, how can I track both boxes in a similar manner?

Thank you very much!

Antoine Rincent, 2019

suffix omp
package omp 1

variable lattconst equal 5.08
variable repetition equal 20
variable temperature1 equal 50
variable temperature2 equal 4200
variable temperature3 equal 1400
variable pression equal 1.0

variable testnum equal 2
variable foldernum equal 9

log EffetNb_{foldernum}_{testnum}.log

---------- Initialize Simulation ---------------------

clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

---------- Create Atoms ---------------------

lattice fcc ${lattconst}

variable demirepetition equal “v_repetition/2”

region whole block 0 {repetition} 0 {repetition} -{repetition} {repetition} units lattice
create_box 1 whole

region box1_1 block 0 {repetition} 0 {repetition} {demirepetition} {repetition} units lattice
region box1_2 block 0 {repetition} 0 {repetition} -{repetition} -{demirepetition} units lattice
region box2 block 0 {repetition} 0 {repetition} -{demirepetition} {demirepetition} units lattice

lattice fcc ${lattconst} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 region box1_1
create_atoms 1 region box1_2
create_atoms 1 region box2

group liquid region box1_1
group liquid region box1_2
group solid region box2

---------- Define Interatomic Potential ---------------------

pair_style eam/alloy
pair_coeff * * Au_GLJ10_3.eam.alloy Au
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

delete_atoms overlap 0.35 all all

velocity all create ${temperature1} 420

---------- Run Minimization ---------------------

variable iterations equal 1000000

variable thermooutput equal “v_iterations/5000”
variable resultsoutput equal “v_iterations/10000”
variable rdfoutput equal “v_iterations/100”
variable dumpoutput equal 400

timestep 0.001
thermo ${thermooutput}
thermo_style custom step temp enthalpy etotal pe ke press vol cella cellb cellc
thermo_modify norm no

min_style cg
minimize 1e-25 1e-25 5000 10000

compute LIQUIDE_liquid liquid rdf 50 * *
compute LIQUIDE_solid solid rdf 50 * *

fix 1 all ave/time {rdfoutput} 1 {rdfoutput} c_LIQUIDE_liquid[] file FixAveTime_liquid_{foldernum}_{testnum}.rdf mode vector
fix 2 all ave/time {rdfoutput} 1 {rdfoutput} c_LIQUIDE_solid[] file FixAveTime_solid_{foldernum}_{testnum}.rdf mode vector

fix 3 liquid npt temp {temperature2} {temperature3} (100*dt) iso {pression} {pression} (1000dt)
fix 4 solid npt temp {temperature1} {temperature3} (100*dt) iso {pression} {pression} (1000dt)

dump mydump all xyz {dumpoutput} EffetNb_{foldernum}_${testnum}.xyz

variable iteration equal “step”
variable tempe equal “temp”
variable enth equal “enthalpy”
variable etot equal “etotal”
variable potential equal “pe”
variable kinetic equal “ke”
variable pressure equal “press”
variable volume equal “vol”
variable cellulea equal “cella”
variable celluleb equal “cellb”
variable cellulec equal “cellc”

fix 5 all print {resultsoutput} "{iteration} {tempe} {enth} {etot} {potential} {kinetic} {pressure} {volume} {cellulea} {celluleb} {cellulec}" file EffetNb_{foldernum}_{testnum}.results screen no title “step temp enthalpy etotal pe ke press vol cella cellb cellc”

run ${iterations}

write_restart *EffetNb{foldernum}_{testnum}.restart

unfix 3
unfix 4

fix 6 all npt temp {temperature3} {temperature3} (100*dt) iso {pression} {pression} (1000*dt)

run ${iterations}

write_restart *EffetNb{foldernum}_{testnum}.restart

Thanks again and have a good day!

2

Hi Antoine,

I am going to comment just on your LAMMPS script and not on your method. It appears that you are setting the velocity of ALL atoms in your system to a distribution representing 50 K at the start. (velocity all create ${temperature1} 420). Then you try to thermostat the liquid part of your simulation to 4200 K? Why not just initialize that part at 4200? Also, is NPT correct here? What have other people done in the literature that have used this approach.

If you want to see separate “temperatures” of your system you must also calculate them. Outputting temp is going to give you the total system temperature; which in your case is going to be misleading. Try outputting the “temperature” of the different groups (consult the documentation for fix ave/chunk and compute temp commands). Another helpful suggestion could be to calculate the ke/atom and visualize your system. All of these steps can help your debug more effectively.

Maybe it is just me, but with systems like these it easier to provide help if you provide an image of your model.
Hope this helps!
Best,
Dylan

3

Hi! The reason why I initialised the temperature to 50K (in this case) to every atoms was because I don’t know what other variables the random script uses other than the seed, so to make sure the temperature is initialised in the same manner everywhere, I put it this way, since the fix NPT temperature ramp has enough iterations for the initially liquid parts, I believe the liquid regions get heated up to 4200K quite fast, and then fall towards the final aimed temperature, whereas the solid part starts at the appropriate temperature and heats up to the final aimed temperature in a more straightforward fashion.

As far as I’ve seen in the literature, most people use a fix NPT to obtain certain temperatures, followed by a fix NVE to verify the system’s state. As this is more of an initial test I’ve skipped this step for simplicity reasons! :slight_smile:

I’ll check out the recommended documentation, I hadn’t seen them and didn’t know how LAMMPS worked to provide the temperature of the system in its usual output, thank you very much!

Sorry for not including an image of the model, considering this will probably end up on the forum where email attachments probably can’t be provided, here’s an imgur link with three “pictures” of the simulation, the first where it begins at the desired “initial state”, where the top and bottom are liquid whilst the middle is solid, the second displays how the solid/liquid interface moves forward towards the middle of the simulated box, the third is once it’s entirely done and the entire simulation box is in a liquid state : https://imgur.com/a/C3Vu5sf

4

here is an alternate idea for how to set up a liquid/solid two phase system.

  • run separate (small) simulations to determine the equilibrium density at the melting point for both solid and liquid phase.
  • set up a simulation box with two areas. it may be worth assigning different atom types (but use the same potential parameters) for the atoms in each part. since you know the equilibrium densities, you can set up the box volume accordingly and don’t need to run a variable cell simulation (which is notoriously tricky and you can properly only do this for the combined system by coupling the one dimension where you have the phase transition)
  • now do simulations to equilibrate each phase separately, i.e. have two groups and only apply time integration to one groups at a time. that will allow you to approach the solid phase from below and without a long time consuming ramp, but starting close to the melting point already. similarly, the liquid phase can be started at a much higher temperature, where melting is easy and the system will rearrange quickly and lose memory of the initial crystal structure. then cool this one down (again, can be rapidly) and equilibrate close to the melting point. you may alternate between time integrating the two parts of the system a couple of times
  • if each part has reached the desired equilibrium, you can switch to applying time integration to the whole system at the expected melting point temperature, while monitoring whether the liquid or solid phase grows and (manually) adjusting the temperature until the phase boundary does not move.

by time integrating each part separately, it is much easier to lead each phase to somewhere near the melting point. the atoms in the other part will keep their velocities, so continuing should be smooth.

axel.

5

Thank you very much for the very well described idea! I’ll look into it with my professor since it seems very promising on a calculation point of view :slight_smile:

Have a good day!
Antoine.