Skip to content
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

Updates to energy minimization routines #1061

Merged
merged 17 commits into from
Dec 5, 2022
Merged

Conversation

chrisiacovella
Copy link
Contributor

@chrisiacovella chrisiacovella commented Nov 12, 2022

PR Summary:

A few extensions to the existing energy minimization functions for compounds. In particular, changes were made to support additional OpenBabel features.

General operations:
1- After energy minimization, shift the Compound to the original COM; this behavior can be toggled with the shift_com variable. Default behavior will be to shift to the original COM. For example, toggle the shift off via:
dodecane = Alkane(n=12) dodecane.energy_minimize(steps=1000,shift_com=False)

2- Rather than retaining the COM position, the user can specify an "anchor" within the Compound. After energy minimization, the Compound is shifted such that the position (or COM as the case may be) of the anchor is unchanged. Example usage:

dodecane = Alkane(n=12) anchor_point = dodecane.labels['chain'].labels['CH3'][0] dodecane.energy_minimize(steps=10000,shift_com=False, anchor=anchor_point)

OpenBabel specific:
1- Energy minimization call will accept a list of Compounds whose positions will be fixed during energy minimization. Note, positions are not truly fixed, but have a strong harmonic energy term associated between the particle and its original position. This can apply to a Compound at any level, as the code will loop over all successors to fix any particle in the Compound. Example usage:

dodecane = Alkane(n=12) point0 = dodecane.labels['chain'].labels['CH3'][0] dodecane.energy_minimize(steps=1000,shift_com=False, fixed_compounds=point0)

Note, can also specify which directions to fix:

dodecane = Alkane(n=12) point0 = dodecane.labels['chain'].labels['CH3'][0] point1 = dodecane.labels['chain'].labels['monomer'][0] fixed = [ [ point0, (True, True, True)], [ point1, (False, False, True)] ] dodecane.energy_minimize(steps=1000,shift_com=False, fixed_compounds=fixed)

2- Distance constraints...a list of Compound pairs and the constrained distance. This must be between Compounds that represent atoms; cannot be between Compounds which contain other Compounds or multiple particles. The code will check to ensure only are applying to Compounds without children. Example usage:

dodecane = Alkane(n=12) point0 = dodecane.labels['chain'].labels['CH3'][0].labels['C'][0] point1 = dodecane.labels['chain'].labels['CH3'][1].labels['C'][0] dodecane.energy_minimize(steps=10000,shift_com=False, distance_constraints=[ [(point0, point1), 0.5 ]])

3- Will accept a list of Compounds to ignore during energy minimization. These Compounds will not have any parameters assigned during energy minimization from the force field, and will have their positions "fixed" (as described above). Note, interactions for ignore particles can be set by a user, e.g., via a distance constraint. This can apply to a Compound at any level, as the code will loop over all successors to fix any particle in the Compound. Example usage:

dodecane = Alkane(n=12) helium = mb.Compound(pos=dodecane.labels['chain'].labels['monomer'][5].pos+[0.0,0.0,0.05], name='He', element='He') dodecane.add(helium) dodecane.energy_minimize(shift_com=False, ignore_compounds=helium)

Can combine ignore_compounds with distance constraints, to, e.g., allow the COM of several atoms to be tied to the position of an associated coarse-grained bead.

PR Checklist


  • Includes appropriate unit test(s)
  • Appropriate docstring(s) are added/updated
  • Code is (approximately) PEP8 compliant
  • Issue(s) raised/addressed?

@codecov
Copy link

codecov bot commented Nov 12, 2022

Codecov Report

Attention: Patch coverage is 90.47619% with 10 lines in your changes missing coverage. Please review.

Project coverage is 90.69%. Comparing base (d051e79) to head (2a2d467).
Report is 141 commits behind head on main.

Files with missing lines Patch % Lines
mbuild/compound.py 90.47% 10 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1061      +/-   ##
==========================================
- Coverage   90.69%   90.69%   -0.01%     
==========================================
  Files          57       57              
  Lines        5397     5500     +103     
==========================================
+ Hits         4895     4988      +93     
- Misses        502      512      +10     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@chrisiacovella
Copy link
Contributor Author

I put together two notebooks demonstrating the usage (these will eventually be merged into the main tutorials). https://github.com/chrisiacovella/mbuild_energy_minimization

Copy link
Contributor

@CalCraven CalCraven left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looking good, such an exciting feature!

mbuild/compound.py Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/tests/test_compound.py Show resolved Hide resolved
mbuild/tests/test_compound.py Show resolved Hide resolved
chrisiacovella and others added 5 commits November 22, 2022 16:51
fixing typo

Co-authored-by: CalCraven <54594941+CalCraven@users.noreply.github.com>
fixing typo

Co-authored-by: CalCraven <54594941+CalCraven@users.noreply.github.com>
fix more typos in comments

Co-authored-by: CalCraven <54594941+CalCraven@users.noreply.github.com>
mbuild/compound.py Outdated Show resolved Hide resolved
Copy link
Member

@daico007 daico007 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, only a few comments to shorten the codes. The logic looks fine, the sanity checks can be a bit disruptive of the flow (may be have those as a helper method?).

mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
if len(fixed_compounds) == 2:
# At this point in the code, we only want to see if the second element in the list
# is a Compound. If it is a Compound we don't need to santize the input.
# If it is not a Compound, (e.g., if the list passed is of the form [Compound, [True, true,True]],
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# If it is not a Compound, (e.g., if the list passed is of the form [Compound, [True, true,True]],
# If it is not a Compound, (e.g., if the list passed is of the form [Compound, [bool, bool, bool]],

mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
Copy link
Contributor

@CalCraven CalCraven left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM! A few minor suggestions I'm seconding from @daico007.

mbuild/compound.py Show resolved Hide resolved
mbuild/compound.py Show resolved Hide resolved
An individual Compound or list of Compounds that will have their
position fixed during energy minimization. Note, positions are fixed
using a restraining potential and thus may change slightly.
Position fixing will apply to all Particles (i.e., atoms) that exist
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are we capitalizing or not capitalizing the word particle? I would maybe suggest no, since we really do want people to understand that these are really Compounds.

mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
mbuild/compound.py Outdated Show resolved Hide resolved
if id(p1) == id(particle):
pid = index
if all_true == True:
ob_constraints.AddAtomConstraint(pid)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this get's fixed by the suggest made hy @daico007 in line 2486 to 2498.

@chrisiacovella
Copy link
Contributor Author

I think I got all the comments and suggestions incorporated.

Copy link
Member

@daico007 daico007 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only one final bit and everything LGTM!

Comment on lines 2488 to 2493
elif dims[0] == True:
ob_constraints.AddAtomXConstraint(pid)
elif dims[1] == True:
ob_constraints.AddAtomYConstraint(pid)
elif dims[2] == True:
ob_constraints.AddAtomZConstraint(pid)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
elif dims[0] == True:
ob_constraints.AddAtomXConstraint(pid)
elif dims[1] == True:
ob_constraints.AddAtomYConstraint(pid)
elif dims[2] == True:
ob_constraints.AddAtomZConstraint(pid)
else:
if dims[0] == True:
ob_constraints.AddAtomXConstraint(pid)
if dims[1] == True:
ob_constraints.AddAtomYConstraint(pid)
if dims[2] == True:
ob_constraints.AddAtomZConstraint(pid)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless they are mutually exclusive, I think they should belong to individual if block

@daico007 daico007 merged commit ee91caa into mosdef-hub:main Dec 5, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants