Dipoles Switching

Hi Everyone,
I wrote this code to get dipoles switching after applying an external electric field, but I’m getting errors on lost atoms after specific time steps. Sometimes at 11 ps and sometimes at 17 ps.
Can you please take a look at my code and tell me what am I doing wrong?
I have changed electric field magnitude and directions, but I’m getting same errors at the same time steps.
My other question is that how can I define specific temperature and keep it fixed during the simulation for sphere atome style?
Best regards,
Hamed

units metal
dimension 3
boundary p p p
atom_style hybrid dipole sphere
atom_modify map array

neighbor 2.0 bin
neigh_modify delay 10 check yes

read_data BiFeO3.dat

pair_style hybrid/overlay coul/streitz 12.0 wolf 0.20 meam/c lj/cut/dipole/cut 10.0
pair_coeff * * meam/c BiFeO3zigzag.library.meam Fe Bi O BiFeO3zigzag.parameter.meam Fe Bi O
pair_coeff * * coul/streitz FeBiO3.streitz Fe Bi O
pair_coeff * * lj/cut/dipole/cut 1.0 1.0

fix 1 all spring/self 1000.0
fix_modify 1 energy yes
fix 3 all qeq/slater 1 12.0 1.0e-6 100 coul/streitz
dump 1 all custom 10000 dump.atom x y z mux muy muz
dump 2 all xyz 10000 dump.xyz
timestep 0.0001
run 10000
fix 5 all efield -0.0003 0.0 -0.0003
fix 6 all nve/sphere update dipole/dlm
compute 1 all property/atom mux muy muz
run 200000
unfix 5
run 200000

Data File:
1 3 1.12238351 1.13820203 1.67962192 -0.81466667 1 1 1 1.2 1.141
2 3 1.12238351 1.13820203 15.7995499 -0.81466667 1 1 1 1.2 1.141
3 3 1.12238351 1.13820203 29.9194779 -0.81466667 1 1 1 1.2 1.141
4 3 1.12238351 1.13820203 44.0394059 -0.81466667 1 1 1 1.2 1.141
5 3 1.12238351 1.13820203 58.1593339 -0.81466667 1 1 1 1.2 1.141

your model/force field looks bogus. where did you take it from?
especially adding an lj/cut term with a giant epsilon of 1eV and a tiny sigma of 1\AA (and the same parameters for all atom types to boot) makes no sense at all.
the fact that you have to use fix spring/self and still are losing atoms, is an indication that your forces are very wrong.

if your model is correct and your simulation settings are proper, your system should conserve energy after it is equilibrated. that is just a basic requirement of good MD simulations. quickly rising temperature is also a sign of a bad force field. your question about a thermostat to force your system to have the desired temperature is a red herring, you are trying to suppress the symptom, but not remove the cause. that will lead to bogus results as well.

GI-GO,
axel.

Hi Axel,
Thanks so much for your reply.

The value that I was using for electrical field was getting from the published paper which it was 1 MV/cm (0.01 V/A - metal units) but in the code that I have sent it to you to avoid lost atoms I was trying smaller value which it is 30 KV/cm(0.0003 V/A), in negative [101] direction.

To integrate the dipole interaction I just used the default command from lammps for bismuth ferrite (BiFeO3) structure. lj/cut/dipole/cut 10.0 * * lj/cut/dipole/cut 1.0 1.0
so are suggesting me to have different paie_coeff for each different pair of my atoms (for example pair_coeff * * lj/cut/dipole/cut 1.0 1.0 and pair_coeff 2 3 1.0 1.0 2.5 4.0) and is there any example or suggestions that I can follow to get the best value for sigmas and epsilons ?

I got your point about temperature so I shouldn’t keep it fixed.

Regards,

You are missing my point.
You are adding meam/c and coul/streitz and lj/cut/dipole/cut on top of each other. that is completely bogus and computing both coulomb and non-coulomb interactions multiple times. adding charge equilibration on top of that is making matters even worse.

i am suggesting that you need to have person-to-person conversation with somebody that understands force fields and how to correctly find and assign such parameters and functional forms to have consistent and correct results. there is no such thing as a “LAMMPS default setting” for force fields. LAMMPS has some examples, but those are just examples.
there is no point at all in your trying to work on interactions of an electric field with dipoles, since those also have to be determined to be consistent with your system.

you are also missing my point about temperature, confirming again, that you are lacking essential training in MD simulations which you need to obtain before continuing with your studies.
there is no simple “do this, not that” answer to set up and run a meaningful simulation of your system, since what you are trying to do is quite advanced and requires an understanding far beyond of what you appear to be trained.

there is no point at all in even trying to study the impact of an electrical field and i have serious doubts whether what you seem to be after can be easily done in LAMMPS without programming work.

axel.