6
$\begingroup$

What I am trying to simulate is an evaporation of ammonia water, after being heated by a gold plate underneath. But I have received a lost atom error during minimization. If I removed the minimization command along with the velocity command, it gives bond atom missing error.

I am assuming that it is due to the water penetrating the gold plate but I have tried out all the fix/wall commands and fix/addforce command, the results are still the same.

Below are my input script:

# Create a simulation box
units real
dimension 3
atom_style      full
bond_style      hybrid morse harmonic
angle_style     harmonic
dihedral_style  charmm
improper_style  umbrella
special_bonds   dreiding
pair_style      hybrid/overlay lj/cut/coul/long 12.0 hbond/dreiding/lj 4 6 6.5 90 
kspace_style    pppm 1.0e-4
kspace_modify   slab 3.0
dielectric      1.0
boundary p p f

region awmix block 0 40.07 0 40.07 -5 100.07  #ammonia water region
create_box 4 awmix bond/types 2 angle/types 2 extra/bond/per/atom 3 extra/angle/per/atom 3 extra/dihedral/per/atom 3 extra/improper/per/atom 3 extra/special/per/atom 3
read_data ATLASparticle.data add append

region plate block 0 40.07 0 40.07 -5 -1 # gold plate
lattice fcc 4.08
create_atoms 4 region plate
lattice none 1.0

group liquid type 2 3 
group solid type 4

# Define interaction potential
pair_coeff 1 1 lj/cut/coul/long 0.0 0.0
pair_coeff 1 2 lj/cut/coul/long 0.00278209 3.0544905
pair_coeff 1 3 lj/cut/coul/long 0.00309354 2.939787
pair_coeff 2 2 lj/cut/coul/long 0.21 3.36 
pair_coeff 2 2 hbond/dreiding/lj 2 i 4.0 2.75 4
pair_coeff 2 2 hbond/dreiding/lj 2 j 4.0 2.75 4
pair_coeff 2 3 lj/cut/coul/long 0.08606498 3.1478565
pair_coeff 2 3 hbond/dreiding/lj 2 i 4.0 2.75 4
pair_coeff 2 3 hbond/dreiding/lj 2 j 4.0 2.75 4
pair_coeff 3 3 lj/cut/coul/long 0.1553 3.166
pair_coeff 3 3 hbond/dreiding/lj 2 i 4.0 2.75 4
pair_coeff 3 3 hbond/dreiding/lj 2 j 4.0 2.75 4
bond_coeff 1 morse 101.905 2.3470 1.0124
bond_coeff 2 harmonic 600.0 1.0
angle_coeff 1 103.045 106.67
angle_coeff 2 75.0 109.47

mass 4 196.96657 # gold

pair_coeff 4 4 lj/cut/coul/long 11.697 2.574

pair_coeff 1 4 lj/cut/coul/long 0.0 0.0
pair_coeff 2 4 lj/cut/coul/long 0.0 0.0
pair_coeff 3 4 lj/cut/coul/long 2.8675 0.0554
pair_coeff 3 4 lj/cut/coul/long 2.8675 0.0114

# Count molecule
variable oxygen atom "type==3"
variable nitrogen atom "type==2"
group oxygen dynamic all var oxygen
group nitrogen dynamic all var nitrogen
variable nO equal count(oxygen)
variable nN equal count(nitrogen)
fix oxy all ave/time 1000 1 1000 v_nO file number_oxygen.dat
fix nit all ave/time 1000 1 1000 v_nN file number_nitrogen.dat

# Count energy
variable kinetic_energy equal ke
variable potential_energy equal pe
variable pressure equal press
variable temperature equal temp
variable volume equal vol
fix myat1 all ave/time 1000 1 1000 v_kinetic_energy v_potential_energy v_pressure v_temperature v_volume file energy.dat

#Provide initial velocity to the atoms by setting temperature of our system = 298.15K
velocity all create 298.15 102939 dist gaussian mom yes rot yes

#Perform minimization
minimize 2.0e-5 0.001 10000 10000
timestep 1
reset_timestep 0

#Setting system
neigh_modify    every 1 delay 0 check yes

# Results
dump dump03 all custom 1000 dump03.lammpstrj id type x y z ix iy iz vx vy vz
thermo_style custom step time temp pe ke etotal press vol
thermo 1000

# run
fix wall solid addforce 1.0 0.0 0.0 region plate
fix 2 all nve/limit 1
fix 3 all temp/berendsen 80 80 1
timestep 0.1
run 50000

I'm not sure if it is related to the data file (ATLASparticle.data = ammonia water data file) but I don't think that's the case because I have run similar simulation (same data file and parameters) and there is no overlapping.

#Create a simulation box
units real
dimension 3
atom_style      full
bond_style      hybrid morse harmonic
angle_style     harmonic
dihedral_style  charmm
improper_style  umbrella
special_bonds   dreiding
pair_style      hybrid/overlay lj/cut/coul/long 12.0 hbond/dreiding/lj 4 6 6.5 90 
kspace_style    pppm 1.0e-4
dielectric      1.0
boundary p p p

region awmix block -50 58.63 -50 58.63 -50 58.63  #ammonia water region
create_box 3 awmix bond/types 2 angle/types 2 extra/bond/per/atom 3 extra/angle/per/atom 3 extra/dihedral/per/atom 3 extra/improper/per/atom 3 extra/special/per/atom 3
read_data ATLASparticle.data add append

group vapor type 2 3

#Define interaction potential
pair_coeff 1 1 lj/cut/coul/long 0.0 0.0
pair_coeff 1 2 lj/cut/coul/long 0.00278209 3.0544905
pair_coeff 1 3 lj/cut/coul/long 0.00309354 2.939787
pair_coeff 2 2 lj/cut/coul/long 0.21 3.36 
pair_coeff 2 2 hbond/dreiding/lj 2 i 4.0 2.75 4
pair_coeff 2 2 hbond/dreiding/lj 2 j 4.0 2.75 4
pair_coeff 2 3 lj/cut/coul/long 0.08606498 3.1478565
pair_coeff 2 3 hbond/dreiding/lj 2 i 4.0 2.75 4
pair_coeff 2 3 hbond/dreiding/lj 2 j 4.0 2.75 4
pair_coeff 3 3 lj/cut/coul/long 0.1553 3.166
pair_coeff 3 3 hbond/dreiding/lj 2 i 4.0 2.75 4
pair_coeff 3 3 hbond/dreiding/lj 2 j 4.0 2.75 4
bond_coeff 1 morse 101.905 2.3470 1.0124
bond_coeff 2 harmonic 600.0 1.0
angle_coeff 1 103.045 106.67
angle_coeff 2 75.0 109.47

#Provide initial velocity to the atoms by setting temperature of our system = 298.15K
velocity all create 298.15 102939 dist gaussian mom yes rot yes

#Perform minimization
minimize 2.0e-5 0.001 10000 10000
timestep 1
reset_timestep 0

#Setting system
neigh_modify    delay 10 every 1 check yes

variable oxygen atom "type==3"
variable nitrogen atom "type==2"
group oxygen dynamic all var oxygen
group nitrogen dynamic all var nitrogen
variable nO equal count(oxygen)
variable nN equal count(nitrogen)
fix oxy all ave/time 1000 1 1000 v_nO file number_oxygen_test.dat
fix nit all ave/time 1000 1 1000 v_nN file number_nitrogen_test.dat

dump dumpMIX all custom 1000 dumpMIXtest.lammpstrj id type x y z ix iy iz vx vy vz
thermo_style custom step time temp pe ke etotal press vol
thermo 1000

variable kinetic_energy equal ke
variable potential_energy equal pe
variable pressure equal press
variable temperature equal temp
variable volume equal vol
fix myat1 all ave/time 1000 1 1000 v_kinetic_energy v_potential_energy v_pressure v_temperature v_volume file energy_test.dat

#run
fix         2 all nve
fix         3 all temp/berendsen 413.15 413.15 100.0
fix         4 all temp/rescale 1 413.15 413.15 0.05 1.0
fix         5 all press/berendsen iso 4.0 4.0 1000.0
timestep 0.1
run 1000000

*You can get the data file from link below: https://github.com/JackyXes/ATLASparticle/blob/main/ATLASparticle.data Please ignore the coordinate values because I changed them for 2 different scripts.

Thanks in advance.

*Edit

enter image description here

As you can see, there are molecules entering the bottom plate. I tried using fix wall/region on the plate region but I got a illegal command error instead.

Improve this question
$\endgroup$
10
  • $\begingroup$ Jack, the link you have posted for the data file doesn't work, maybe because it's a private repository? Please add the data file with your question, so that I can reproduce your issue. $\endgroup$
    – Magic_Number
    Nov 1, 2022 at 6:14
  • $\begingroup$ Also, an observation from the codes, in the 1st case you have used fixed boundary conditions in the z-direction. If an atom moves outside the box from +/- z direction it will be deleted on the next timestep (unlike in periodic boundary conditions where it is reintroduced again from the other direction). This will typically generate an error unless you have set the thermo_modify lost option to allow for lost atoms. $\endgroup$
    – Magic_Number
    Nov 1, 2022 at 6:17
  • 1
    $\begingroup$ HI @Magic_Number, after running with more recorded timestep, I think the main reason is because the molecule pass through zlo and have atom deleted, as result, the bond connecting that atom is lost. But I know having a fixed z boundary is possible to run because that's what I found from several lammps papers. Just that I am missing a command to prevent the molecules from passing through the zlo boundary. Once again, thanks. $\endgroup$
    – Jack
    Nov 1, 2022 at 8:32
  • 1
    $\begingroup$ Let us continue this discussion in chat. $\endgroup$
    – Magic_Number
    Nov 1, 2022 at 9:34
  • 3
    $\begingroup$ @Magic_Number next time please use the LAMMPS chat room instead of creating a new one. Please support this proposal if you feel that this process should be simpler. $\endgroup$
    – Nike Dattani
    Nov 1, 2022 at 16:23

1 Answer 1

Reset to default
6
$\begingroup$

You have made many changes between the second attached input (which is a bulk ammonia-water mixture at 413K -- albeit with a very strange choice of two thermostats) and the first input. The first input has an ammonia-water mixture:

  • with one non-periodic axis
  • in an extremely inhomogeneous box (there's vacuum!)
  • and a slab of gold
  • all thermostatted together
  • at 80K

Any of those changes, or combination of those changes could have broken your simulation (and remember -- just because your first simulation ran does not even guarantee that it is correct). You have to test them one by one and proceed via elimination.

(Note that "my input crashed" questions are usually not easily answered on StackExchange because the diagnostic process invariably involves some back-and-forth. The process is more suited to a forum, such as the LAMMPS forum at MatSci.org.)

Improve this answer
$\endgroup$
11
  • $\begingroup$ Actually I posted the second script is just to show that the data file is not the one causing this issue. What I trying to do is to simply prevent the molecules from passing through the gold plate. $\endgroup$
    – Jack
    Nov 1, 2022 at 10:03
  • $\begingroup$ Saying what does not cause the issue is of limited help in working out what does cause the issue. You have one "working" input which, after a few changes, yields one "non-working" input -- I simply observe that almost certainly one of those changes was the one that made your input not work, and finding that one change is what you might want to do. $\endgroup$
    – Shern Ren Tee
    Nov 1, 2022 at 10:20
  • $\begingroup$ And if what you are trying to do is to "simply prevent the molecules from passing through the gold plate", then how do you know if you're still making a useful model of a physical system, not just filling a hard drive with numbers? After all, you could charge all your particles negative and let them repel each other, or accelerate all your molecules upwards at 2 g, or fill a plexiglass box with ball bearings. (I'm not being rude. This is an actual experiment we ran as second year undergrads to observe the barometric distribution.) $\endgroup$
    – Shern Ren Tee
    Nov 1, 2022 at 10:26
  • 2
    $\begingroup$ If you comment too much the system will try to move you to chat. If that happens please don't click the button. We don't want more chat rooms right now. Please use an existing one. $\endgroup$
    – Nike Dattani
    Nov 1, 2022 at 13:24
  • 1
    $\begingroup$ Stack Exchange tends to handle these kinds of questions by back-and-forth in the comments and then posting the eventual solution as an answer. $\endgroup$
    – user253751
    Nov 1, 2022 at 14:58

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .