-
Notifications
You must be signed in to change notification settings - Fork 26
Added simple rounding error solution #1167
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Changes from 2 commits
dda5d58
2256fb8
e9bfd28
d6cbe74
84922e3
9fb2794
a8223e6
dd48e2e
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -222,7 +222,38 @@ def _write_atoms( | |
| top.write("[ atoms ]\n") | ||
| top.write(";index, atom type, resnum, resname, name, cgnr, charge, mass\n") | ||
|
|
||
| for atom in molecule_type.atoms: | ||
| def _adjust_charges(charges: numpy.array, tolerance=8) -> numpy.array: | ||
| """ | ||
| Adjust charges so that written charge for a molecule type is 0. | ||
| """ | ||
| rounded_charges = numpy.round(charges, tolerance) | ||
| total_charge = numpy.round(numpy.sum(charges_to_write), 0) | ||
|
|
||
| def _rounding_error(_arr, _sum, _tolerance): | ||
| return numpy.round(numpy.sum(_arr) - _sum, _tolerance) | ||
|
|
||
| rounding_error = _rounding_error(rounded_charges, total_charge, tolerance) | ||
|
|
||
| if rounding_error != 0: | ||
| rounded_charges += numpy.round( | ||
| -rounding_error / len(charges_to_write), | ||
| tolerance, | ||
| ) | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. (blocking) Wouldn't 0 charges be affected here, contradicting the "never adjust 0 charges" comment below?
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. You are right, should be resolved now. |
||
| diff = _rounding_error(rounded_charges, total_charge, tolerance) | ||
| while diff != 0: | ||
|
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. (not blocking) Under what circumstances would this while loop iterate more than once? I can't think of any.
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I agree, since tolerance is exactly the same, then one round should be enough. |
||
| # Never adjust 0 charges | ||
| indices = numpy.where(rounded_charges != 0)[0] | ||
| rounded_charges[indices[0]] = numpy.round( | ||
|
pbuslaev marked this conversation as resolved.
Outdated
|
||
| rounded_charges[indices[0]] - diff, | ||
| tolerance, | ||
| ) | ||
| diff = _rounding_error(rounded_charges, total_charge, tolerance) | ||
| return rounded_charges | ||
|
|
||
| charges_to_write = numpy.array([atom.charge.m for atom in molecule_type.atoms]) | ||
| rounded_charges = _adjust_charges(charges_to_write) | ||
|
|
||
| for atom, charge in zip(molecule_type.atoms, rounded_charges): | ||
| if merge_atom_types: | ||
| top.write( | ||
| f"{atom.index:6d} " | ||
|
|
@@ -231,7 +262,7 @@ def _write_atoms( | |
| f"{atom.residue_name:8s} " | ||
| f"{atom.name:6s}" | ||
| f"{atom.charge_group_number:6d}" | ||
| f"{atom.charge.m:20.12f}" | ||
| f"{charge:20.12f}" | ||
| f"{atom.mass.m:20.12f}\n", | ||
| ) | ||
| else: | ||
|
|
@@ -242,7 +273,7 @@ def _write_atoms( | |
| f"{atom.residue_name:8s} " | ||
| f"{atom.name:6s}" | ||
| f"{atom.charge_group_number:6d}" | ||
| f"{atom.charge.m:20.12f}" | ||
| f"{charge:20.12f}" | ||
| f"{atom.mass.m:20.12f}\n", | ||
| ) | ||
|
|
||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
(not blocking) My external perspective here is that needing to come back to this definition repeatedly as this is called in multiple places decreases readability substantially, and that the line length is about the same if
numpy.roundwere called directly on each line.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I see the point, but when initially I coded it without a function the code looked to heavy as for my liking. With the function it looked tidier, so I opted for having this function. I don't have any strong opinion here, so I am happy to get rid of the function if this is more in line with openff codestyle.