Closed lazaroid closed 4 years ago
Thanks for that. As you see, in that line, I've commented out the round value you've suggested. Before reverting it, I need to re-investigate why I did this.
I suspect it was simply a case of combining the relevant fundamental constants to get your own, more precise value. Amusingly, the value 18.22223 (one extra 2...) crops up in a couple of comments within AmberTools in a couple of places!
Typically, I create gromacs topologies using a combination of antechamber, tleap, and then acpype.py for the final step. In generaly, the topology will be determined for a single, small (<100 atom), uncharged molecule before I combine a few hundred of them into a simulation box for a condensed-phase simulation. Occasionally, a small amount of charge creeps in and gromacs complains.
I eventually tracked this down to the value of qConv. I forget whether it is antechamber or tleap (I originally worked this out several months back and submitted it to what I thought was the "official" github repository for acpype, before the owner pointed me at this one, after deleting my issue!) but atomic charges are converted from units of "e" to kCal, and then acpype.py converts these charges back into "e" units.
Within AmberTools, the conversion factor is exactly 18.2223; in acpype.py, the conversion factor used is the more precise:
qConv = 18.222281775 # 18.2223
(line 129).More precision would normally be "better" but in this case, conversion in one direction is less precise than in the other, and this difference can sometimes give a change in the 5th or 6th decimal place (again, I forget which) that is too small for the total charge cutoff value checked for by gromacs (1e-4). However, once a few hundred of these molecules are added to the system, the cutoff value is exceeded and gromacs gives a warning about the system having an overall charge.
Changing qConv to be 18.2223 completely removes this issue.