Dear Lammps,
I am simulating a system of rigid linear molecules. When simulated in the bulk (cubic volume with periodic boundary conditions) everything works fine. When I introduce an interface by changing zlo,zhi in the read_data file, I find that the fix rigid integrator breaks down and for some molecules, the molecule end-to-end distance D is no longer conserved. Thus these molecules do not behave like rigid molecules. The fix rigid integrator fails to conserve energy and my kinetic energy blows up. Potential energy does not behave too strangely.
I've done extensive testing on the bulk and interfacial systems. I find that the molecules whose D values are not conserved started out lying across one of the boundaries, either x or y. For the z-boundary, since I introduce an interface there by changing zlo, zhi, I go through the initial config file and adjust the positions so that the molecule is not cut off by the boundary. I do not do this for the x and y boundaries since these are periodic. For the bulk system (cubic volume, no interface), D is conserved for all molecules within a window of +/- 1% of D. Upon introducing an interface to the same config, correcting any molecules lying across a z-boundary, I find that moleculels have D values outside of D+/- 1%, and that these systems do not conserve energy.
Rigid body molecule dynamics should conserve D within a window much smaller than +/- 1%. I think either I am violating some implicit lammps assumption about fix rigid or there is a bug with this integrator when interfaces are introduced.
Any comments and suggestions would be greatly appreciated!
Thanks,
Joyce