I have a confusion about the build-up of the simulation model for lammps. Currently, I am interested in simulating the dipole moment of centrosymmetric crystals (i.e. MOFs with hundreds or thousands of atoms, bonds etc.). Theoretically, due to the centrosymmetry, the total dipole moment of the entire model should be zero. however, I found a shortcoming of the model building, that is there is no program or software or even lammps commands to restrict the structure to keep its centrosymmetry during simulation, thus the resulting total dipole moment is not zero.
I have tried to use the programs msi2lmps.exe and cif2lammps to convert the file to lammps format, but both of them cannot handle the symmetry operation, so I removed the symmetry to make the model in P1 space group before conversion. and then, the orientation and magnitude of the dipole moment of two centrosymmetric parts (atoms) can not eliminate each other.
does anyone know if there is a way to restrict the model to keep the centrosymmetry during the simulation? and please correct me if I made the wrong description. thank you so much.
That would be only possible at 0K temperature. At finite temperature, even in a perfect crystal, atoms would vibrate around their crystal positions and thus the total dipole would not be exactly zero.
If you start from a perfectly symmetric geometry and do not have any operations that break this symmetry, it will be preserved since the symmetric geometry will result in symmetric forces. If you assign velocities from a random distribution or due to floating point math limitations (addition in floating point is not associative so the order of summing can lead to tiny differences in the total result) you can slowly diverge from the perfectly symmetric behavior. The same will happen in real situations due to thermal noise. Now it depends on the intrinsic behavior of the model whether it is a stable configuration and will preserve symmetry and pull diverging atoms back, or whether it is a metastable configuration and the atoms will relax into a more stable configuration.
LAMMPS will just follow the forces and use the provided geometry, atom typing, atom charges and so on. So whether the symmetry is retained is a matter of using a suitable force field, applying it correctly, and starting from a geometry that is consistent with the force field. Any additional “guiding force” would be entirely artificial and thus your simulation would turn from a scientific experiment into a computer animation. Thus it would not be a solving a problem but just suppress it.
As already explained, the total dipole moment cannot be exactly zero.
Thank you so much, Axel, your kind answer helps a lot.
I am trying to understand the way LAMMPS calculates the dipole moment. there is a direct way to obtain a dipole moment - “compute dipole” command, and I can use fix ave/time to obtain .
but I also tried to use another way to calculate it:
M=sum (qi*ri)
here I only show the dipole moment along the x direction, the commands I used are:
compute peratom all property/atom xu q
variable Mperatom atom “c_peratom[1]*c_peratom[2]”
compute sum all reduce sum v_Mperatom
fix lMxl all ave/time {nevery} {nrepeat} ${nfreq} c_sum
the problem is the values obtained in two ways are quite different, can you give me any suggestions to revise it?
No. What compute dipole does is explained in the documentation.
If you need to know more details, you need to look at the source code.
I don’t know your system and how different those numbers are.
If you want to debug your LAMMPS computation, you need to set up some example test systems to validate that your computation is correct. I have no time to do that for you.
sorry for providing the missing info. My system is cubic and centrosymmetric (MOF: HKUST-1), so it is isotropic along each direction, and the data file contains bonds, dihedral and other interactions. By using “compute dipole”, the obtained value is zero, while for the second way, the value is around -108, this negative value also confused me.