Chain of particles starts to rapidly rotate

[NicholasHattrup]

Hello all,

I am trying to measure the thermal conductivity along the x-direction for a group of 3 ‘nano-particles’ of Bi2Te3 with periodic boundaries along x & y and slab boundary along z. The workflow for this system set up is as follows:

1. Generate particle chain and relax x-direction via NPT with a fix spring tether command for each nanoparticle group to keep their y and z com coordinates the same throughout relaxation along the x-axis
2. Relax in NVT while removing the tether to ensure chain slowly relaxes followed by equilibration in NVE before production phase
3. Production NVE where I measure the heat-flux and then compute thermal conductivity via the green-kubo relations

Below are my NPT and NVT relaxation scripts:

NPT:

units metal 
dimension 3 
boundary p p f 
atom_style charge
processors 3 2 1
pair_style hybrid/overlay morse 5.5 coul/long 12 


neigh_modify every 1 delay 0 check yes page 100000 one 10000
neigh_modify exclude none


pair_modify shift yes
kspace_style pppm 1.0e-4
kspace_modify slab 3.0

read_data chain_3.data

# type 1 = Bi
# type 2 = Te1
# type 3 = Te2


pair_coeff * * coul/long 
pair_coeff 1 1 morse 0.085 2.212 4.203 5.5 # type1 type2 D a r_0 r_cutoff
pair_coeff 2 2 morse 0.076 1.675 3.642 5.0
pair_coeff 3 3 morse 0.066 2.876 4.312 5.0
pair_coeff 1 2 morse 0.975 1.285 3.089 4.0
pair_coeff 1 3 morse 0.582 1.257 3.251 4.0
pair_coeff 2 3 morse 0.807 0.731 4.497 5.5

timestep 0.0005 # 0.5 femtoseconds 
thermo 5000

variable ke equal ke
variable pe equal pe
variable etotal equal etotal
variable volume equal vol
variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz
variable pxx equal pxx
variable pyy equal pyy
variable pzz equal pzz



# Groups for particle chains

group 1 id 1:473
group 2 id 474:946
group 3 id 947:1419

fix pull_1 1 spring tether 1.0 NULL 22.3457853 31.32114068 0
fix pull_2 2 spring tether 1.0 NULL 22.3457853 31.32114068 0
fix pull_3 3 spring tether 1.0 NULL 22.3457853 31.32114068 0


fix track_energies all ave/time 1 1 20000 v_ke v_pe v_etotal file tether_npt_energies.dat

dump mydump all custom 200000 tether_npt_configs.dat type q x y z # every 100 ps 
dump_modify mydump sort id 




fix my_npt all npt temp 300.0 300.0 $(500.0*dt) x 0 0 $(5000*dt) # Slower equilibration

# Zero net linear and angular momentum 
velocity 1 zero linear
velocity 1 zero angular
velocity 2 zero linear
velocity 2 zero angular
velocity 3 zero linear
velocity 3 zero angular
compute my_stress all pressure thermo_temp


fix track_stress all ave/time 400 10 20000 c_my_stress[1] c_my_stress[2] c_my_stress[3] c_my_stress[4] c_my_stress[5] c_my_stress[6] file tether_npt_stresses.dat

run 10000000 # 5 ns
  
write_data rNEMD_npt_equil.data

NVT:

units metal

dimension 3

boundary p p f

atom_style charge

pair_style hybrid/overlay morse 5.5 coul/long 12

neigh_modify every 1 delay 0 check yes page 100000 one 10000

neigh_modify exclude none

pair_modify shift yes

kspace_style pppm 1.0e-4

kspace_modify slab 3.0

read_data rNEMD_npt_equil.data

# type 1 = Bi

# type 2 = Te1

# type 3 = Te2

pair_coeff * * coul/long

pair_coeff 1 1 morse 0.085 2.212 4.203 5.5 # type1 type2 D a r_0 r_cutoff

pair_coeff 2 2 morse 0.076 1.675 3.642 5.0

pair_coeff 3 3 morse 0.066 2.876 4.312 5.0

pair_coeff 1 2 morse 0.975 1.285 3.089 4.0

pair_coeff 1 3 morse 0.582 1.257 3.251 4.0

pair_coeff 2 3 morse 0.807 0.731 4.497 5.5

timestep 0.0005 # 0.5 femtoseconds

thermo 5000

variable ke equal ke

variable pe equal pe

variable etotal equal etotal

variable volume equal vol

variable pxy equal pxy

variable pxz equal pxz

variable pyz equal pyz

variable pxx equal pxx

variable pyy equal pyy

variable pzz equal pzz

# Groups for particle chains

group 1 id 1:473

group 2 id 474:946

group 3 id 947:1419

fix pull_1 1 spring tether 1.0 NULL 22.3457853 31.32114068 0

fix pull_2 2 spring tether 1.0 NULL 22.3457853 31.32114068 0

fix pull_3 3 spring tether 1.0 NULL 22.3457853 31.32114068 0

fix track_energies all ave/time 1 1 20000 v_ke v_pe v_etotal file tether_nvt_energies_1.dat

dump mydump all custom 200000 tether_nvt_configs_1.data type q x y z

dump_modify mydump sort id

fix my_nvt all nvt temp 300.0 300.0 $(1000.0*dt) # Slower equilibration

# Zero net linear and angular momentum

velocity 1 zero linear

velocity 1 zero angular

velocity 2 zero linear

velocity 2 zero angular

velocity 3 zero linear

velocity 3 zero angular

run 200000 # 100 ps

unfix pull_1

unfix pull_2

unfix pull_3

velocity all create 300.0 704812 # Regenerate velocities

# Zero net linear and angular momentum

velocity 1 zero linear

velocity 1 zero angular

velocity 2 zero linear

velocity 2 zero angular

velocity 3 zero linear

velocity 3 zero angular

variable ang_1 equal angmom(1,x)

variable ang_2 equal angmom(2,x)

variable ang_3 equal angmom(3,x)

compute my_stress all pressure thermo_temp

fix track_angmoms all ave/time 2000 10 20000 v_ang_1 v_ang_2 v_ang_3 file tether_angmoms_1.dat

fix track_stress all ave/time 2000 10 20000 c_my_stress[1] c_my_stress[2] c_my_stress[3] c_my_stress[4] c_my_stress[5] c_my_stress[6] file tether_stresses_1.dat

run 12000000 # 6 ns

write_data rNEMD_nvt_equil_1.data

What is most curious about my set-up is that it seems depending on what velocity seed I use, I can generate initial conditions that lead to the chain rapidly rotating about the x-axis, with increasing rotations over the trajectory.

Here is the angular momentum plotted throughout the NVT trajectory over 12 random-seeds, with the average across the 12 seeds in dashed-red, each separate plot corresponds to one of the three particles in the 3-particle chain:

Which lead to widely different final thermal-conductivity estimates:

The 0 estimates along y and z give me confidence my simulations are not fundamentally incorrect as the finite length along both directions for the chain necessarily must give 0 thermal conductivity, so I have no clue why 2 of the 12 random seeds, with one in particular, absolutely blows up due to whatever random velocity seed I assigned to it.

Has anyone ever experienced a problem like this before? My first thought was there must be some type of in-phase perturbation that continues to grow throughout the simulation time, but I wasn’t able to find anything similar, albeit I don’t really know the correct terminology to describe what I’m seeing. Thanks for any insight someone is able to provide on this!

[akohlmey]

@NicholasHattrup

“”"

FYI, the marker for quoted text is three backquotes not three double quotes.
This is not a python multi-line string, but a markdown syntax element.

Does your input text look like it is properly quoted to you?

[NicholasHattrup]

Thanks for the heads up. I’ll be sure to format properly in the future. Are you recommending to resubmit with the edited formatting for better legibility?

[akohlmey]

You can just edit the existing post

[Germain]

Hi @NicholasHattrup,

@NicholasHattrup

velocity 1 zero linear
velocity 1 zero angular
velocity 2 zero linear
velocity 2 zero angular
velocity 3 zero linear
velocity 3 zero angular

Is there any reason to proceed using these groups? Since all your particle have the same dynamics in each simulation, I suspect a common remnant angular momentum to be left in the system. There seems to be remnant angular momentum on your particle at the start of your data acquisition.

Have you tried adding velocity all zero angular using the default all group?

[NicholasHattrup]

I have not, I’ll try that now on the seed that blew up, but would you be able to point me to anything I can look into to understand this a bit better?

In my mind, if I am zeroing the 3-particles separately, then at-least along the same axis, it should be zero within computer precision as if one had zero-ed out all particles. This clearly would work for linear momentum, so I’m struggling to see the difference for angular. Again I note it depends on the axis one is rotating through.

And even if this didn’t zero it out at the start, is there any reasoning for why there is such a dramatic drift upwards in some of the seeds depending on what velocity seed one starts with?

Thanks so much!

[Germain]

It is hard to give meaningful help and relevant comments without any input file or minimal working example reproducing the behavior.

@NicholasHattrup

And even if this didn’t zero it out at the start, is there any reasoning for why there is such a dramatic drift upwards in some of the seeds depending on what velocity seed one starts with?

My 2cts: there can be several reasons for that. One can be bad luck and ending up with some resonance with your NVT NH thermostat which ends up dramatically increasing a remnant angular momentum. So another idea could be to modify your thermostat and see if you see any correction. Either modifying the tchain parameter or switching to a langevin+nve integrator. You could also start by checking if you can see such effects with a single nve integrator.

[srtee]

Without any bonds, how are your nanoparticles being held together? If your force field somehow supports a stable lattice structure, wouldn’t three nanoparticles be unstable against coalescing into a single spherical nanoparticle?

[NicholasHattrup]

My trajectories are at most 10 nanoseconds, so it very well may be that my structures are metastable and would coalesce over longer time-scales. The potential was pulled from this paper:
https://doi.org/10.1103/PhysRevB.80.165203 . They mention they fit the potential using the “total-energy surface” obtained via ab initio calculations. From here they used this potential to calculate the optimal lattice structure and measure bulk thermal conductivity.

My understanding is that it was validated for the lattice but not strictly fit to only match that type of configuration for the material, one of those try on a nanoparticle and see how it performs as opposed to generating a whole new potential … so far no glaring issues with the potential itself and enough of a working set up to proceed forward with thermal-conductivity calculations.

[srtee]

Have you visualised your trajectories?

[NicholasHattrup]

Yep! Particles remain defined throughout the runs. There is a bit of interfacial ‘smooshing’ between particles, but for now that’s fine. Here is what the 3 particles look like at the start of the run, and about ~11 ns into equilibration. I realize 11 ns is very long, but I am trying to get the residual oscillations damped out prior to proceeding to a second equilibration in NVT and final NVE production run.

I also monitor the stresses throughout these runs as well, I equilibrated over two different trajectories via different seeds for initial velocities (same starting configuration). Below are what those results look like so far along with damping/drag parameters in each “phase” of the NPT equilibration. Blue is the stress in the x-direction, red is the length of the box in the x-direction:

[srtee]

Well, the nanoparticles certainly remain distinct even as they touch. I don’t know MEAM MD well enough to say meaningfully whether the simulation appears valid or not beyond that.

I will say that since periodic boundary conditions result in a system without global rotational symmetry; therefore, I’m not sure there is any guarantee that the total angular momentum is conserved. (Consider an infinitely periodic box with one electric dipole whose moment has constant magnitude. The configuration has some non-zero self-energy due to the dipole’s interactions with its images. Translating the dipole will leave this energy constant. Rotating the dipole will, in general, not leave this energy constant.)