Hello lammps people,
I have been working with lammps for the past month or so and this is my first email to the list.
I have implemented some simple potentials to simulate a very coarse dna model. I have implemented a dihedral, a hard-core repulsive pair interaction and a harmonic-cosine breakable bond potential.
My simulations have ran on single processors for many days with no problems and produced consistent results. Now I am attempting to do parallel simulations with no luck. Any help is appreciated. Here is as much info as I believe is necessary.
I am using 30 Oct stable version + my 3 potentials.
My simulations run fine on single processors for 10+ days.
Using real units
There are 1600 particles in this particular system
mpirun -np 6 lmp_murat -in mysim.in
ERROR on proc 3: Bond atoms 560 1041 missing on proc 3 at step 0 (../neigh_bond.cpp:65)
mpirun -np 4 lmp_murat -in mysim.in
ERROR on proc 0: Bond atoms 301 1300 missing on proc 0 at step 155614 (../neigh_bond.cpp:65)
lmp_murat -in mysim.in > runs fine!
So the error happens in first step for 6 cpus, in some arbitrary step for less cpus. In any case there are more than 100 particles per cpu.
What I have tried and did not work:
Dialing the timestep wayy down : from 75 which works nicely on single processor to 5
frequent neighbor list rebuilds : neigh_modify delay 0 every 1 check no
large neighborlist cutoff : neighbor nsq 50 (from 5 which works on single processor)
alternative decompositon : balance 1.2 rcb with comm_style tiled
3rd law : newton on/off
The catch :
The relevant (missing) atoms are bound by harmonic-cosine bonds which I implemented myself.
My harmonic cosine bond is breakable (it has a cutoff distance beyond which it does not do anything), which means bonded atoms can be arbitrarily distant. This bond is in fact not only breakable but also re-bondable which means I cannot just remove the bonds once they break. This is not a problem in a single processor run, but I am beginning to think that it could be for parallel runs.
I would appreciate any insight into the problem and suggestions to better implement such a breakable bond. I am attaching the source files for my 'hacks'.
Thank you very much
bond_harcos.cpp (6.68 KB)
dihedral_breakable.cpp (11.2 KB)
pair_hardcore.cpp (9.18 KB)