gmx density hangs when computing charge density with -symm
gmx density hangs when computing charge density in the system with shell particles if -symm is provided:
$ gmx density -dens charge -f traj_comp.xtc -sl 100 -symm ... Select the group to center density profiles around: Group 0 ( System) has 2292 elements Group 1 ( Other) has 225 elements Group 2 ( AU) has 225 elements Group 3 ( Water) has 2067 elements Group 4 ( SOL) has 2067 elements Group 5 ( non-Water) has 225 elements Select a group: 2 Selected 2: 'AU' Select 1 group to calculate density for: Group 0 ( System) has 2292 elements Group 1 ( Other) has 225 elements Group 2 ( AU) has 225 elements Group 3 ( Water) has 2067 elements Group 4 ( SOL) has 2067 elements Group 5 ( non-Water) has 225 elements Select a group: 4 Selected 4: 'SOL' Reading frame 0 time 0.000 ...HANGS FOREVER...
If no -symm is asked it works fine.
The system has shell particles for polarization, which may cause the bug. In other systems without shells this option works well.
#1 Updated by klark chen about 3 years ago
My gmx density (gmx 5.1.4) also hangs when computing charge density on a bilayer system with option -center. Since -symm automatically turns on -center(as stated in the document), it is most probably a reproduction of this bug, and codes related to -center could be traced in priority.
#2 Updated by klark chen about 3 years ago
The issue relevant to -center and charge density was mentioned by Reid Van Lehn in Feature #1168. When charge or number is provided to -dens, the code simply overrides atom mass of the topology by atom charges or 1 before COM is calculated. Consequently center of charge or number is actually used for centering the system, which may not be what users desire. On the other hand, in typical systems the total charge can be very close to zero, and center_coords, the function used to calculate the center is designed for center of mass, which applies divisions by the total mass/charge, thus numerical instability may be introduced.
I made a fix for gmx_density.c to center the system by COM in all -dens cases. Quantities for density calculation are hold in a new array den_val to avoid direct overwriting of atom mass.