Dear LAMMPS users and dear Steve,
Thanks for considering my question.
I was finally able, with some external help, to solve the problem!
The reason for this behavior was the data file: indeed, by looking at the z coordinates of the particles in “data”, one can see that some atoms have z coordinates very close to zlo=0, and other atoms have z coordinate very close to zhi=0.1 (I am using a finite but narrow zlo-zhi range, as suggested on the manual https://lammps.sandia.gov/doc/Howto_2d.html). The latter particles also have a z flag value of -1.
It was surprising to me to discover this problem with the data file, and I’ll explain why below (although probably I should do so in a separate post…)
When I start the simulation, I assign to each atom a z coordinate z=0.05=(zhi-zlo)/2 (again, as suggested on the manual: “For each atom in the file, assign a z coordinate so it falls inside the z-boundaries of the box”).
The initial coordinates of the molecules (pentagons) are random, so that initially overlap are present.
In order to remove these overlaps, I use the script below, in which two consecutive runs are performed with a soft potential (pair_style soft).
The strength of the potential (“A” constant) is ramped from 0 to 20 during the first (shorter) run, then from 20 to 100 during the second (longer) run.
The data file written at the end of this simulation (through the write_data command) is then used as input for the script reported in the original post.
When running this script on a data file in which all the atoms have z=0.05 (see attached file “data_initial”), the result is the change of value of the z coordinates described above.
However, it is sufficient to set z=0 for all particles for this problem to disappear.
If this is done, the energy drift problem described in my initial email also disappears!
It is not clear to me what causes the script below to give different results depending on whether the atoms initially have z=0.05 or z=0.
I attach two data files, one in which all the atoms have z=0 and one in which they have z=0.05. In both cases, zlo=0 and zhi=0.1.
Below is the script used to remove the overlaps.
You can verify that the data file written at the end of the simulation is different depending on which file is used as input.
Thanks again for your help.
Best regards,
Valerio
PS: Another interesting (?) observation: when simulating the pentagons composed by 5 atoms (in which no intra-molecular overlaps are present – file data_5_initial attached below), the z coordinate of all the atoms is set to 0 after the script to remove the overlap (below) is run (!). This behavior is different from the one found for the pentagons composed by 10 atoms, which I described above.
I meant to dig a bit more into this problem, to see whether it was possible to reproduce the same behaviors with simpler shapes instead of pentagons, but I still haven’t had the time…
SCRIPT USED TO REMOVE OVERLAPS:
variable tpushoff_1 equal 1e4
variable tpushoff_2 equal 1e6
variable mytemp equal 1.0
variable thermodump equal ((v_tpushoff_1+v_tpushoff_2)/1e3)
variable trajdump equal ((v_tpushoff_1+v_tpushoff_2)/20)
variable tstep equal 0.0001
variable damp equal 1.0
variable soft_cutoff equal 1.1225
variable soft_energy_max_1 equal 20.0
variable soft_energy_max_2 equal 100.0
variable lseed equal 34314142
variable fname string data
variable thermofile string thermo.dat
variable simname string push
units lj
boundary p p p
atom_style bond
dimension 2
read_data ${fname}
pair_style soft {soft_cutoff}
pair_coeff * * 0.0 {soft_cutoff}
variable soft_energy equal ramp(0,${soft_energy_max_1})
neighbor 0.5 bin
neigh_modify every 5 delay 0 check yes
neigh_modify exclude molecule/intra all
variable step equal step
variable temp equal temp
variable etot equal etotal
variable ke equal ke
variable pe equal pe
variable press equal press
thermo_style custom step etotal ke pe temp press
thermo ${thermodump}
timestep ${tstep}
fix fix_print all print {thermodump} "{step} {etot} {ke} {pe} {temp} {press}" file {thermofile} screen no title “#1:step 2:etot 3:ke 4:pe 5:temp 6:press”
fix fix_rigid_damp all rigid molecule langevin {mytemp} {mytemp} {damp} {lseed}
fix fix_zero_linmom all momentum 10 linear 1 1 1
fix fix_zero_angmom all momentum 10 angular
fix fix_push_1 all adapt 1 pair soft a * * v_soft_energy
fix fix_2d_1 all enforce2d
run ${tpushoff_1}
unfix fix_push_1
unfix fix_2d_1
variable soft_energy equal ramp({soft_energy_max_1},{soft_energy_max_2})
fix fix_push_2 all adapt 1 pair soft a * * v_soft_energy
fix fix_2d_2 all enforce2d
run ${tpushoff_2}
write_data data.{simname}
write_restart restart/restart.{simname}
data_initial_z0 (12.4 KB)
data_initial (12.4 KB)
data_5_initial (6.21 KB)