-
Notifications
You must be signed in to change notification settings - Fork 202
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
T019 - Molecular dynamics simulation #56
Conversation
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
Hey @jaimergp, I could not render with sphinx easily because this notebook is probably missing in some sort of configuration file. How can I include it? |
Check docs/talktorials and create the correaponding nblink file there! |
@AndreaVolkamer Everything works well as notebook, colab notebook or html 🎉 Ready for review! |
@Mika-Le - your talktorial's final index will be |
@@ -0,0 +1,1155 @@ | |||
{ |
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.
@@ -0,0 +1,1155 @@ | |||
{ |
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.
Prepare the protein ligand comples -> complex
Suggestion: I would move the 'Merge protein and ligand' one level lower:
- Prepare the protein ligand complex
- Prepare protein
- Prepare ligand
- Merge protein and ligand
Reply via ReviewNB
@@ -0,0 +1,1155 @@ | |||
{ |
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.
You so nicely explained what to expect in all references, despite
Pierre-Simon Laplace, Oeuvres Complètes de Laplace. Théorie Analytique des Probabilités (volume VII Gauthier-Villars (1820), 3rd ed)
Could you add sth here also?
Reply via ReviewNB
@@ -0,0 +1,1155 @@ | |||
{ |
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.
@@ -0,0 +1,1155 @@ | |||
{ |
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.
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 thought I remembered from our first Hackathon that we want to avoid HTML Syntax, which is the only way I know to scale images in notebooks. But I cannot find the guidelines anymore and might be mistaken. Maybe Jaime can tell us.
@@ -0,0 +1,1155 @@ | |||
{ |
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.
@@ -0,0 +1,1155 @@ | |||
{ |
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.
This is starting to go into the details (and maybe too much for an introductory talktorial), but yould we at least link to where to find more information about what e.g. the 'LangevinIntegrator' does?
Reply via ReviewNB
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.
Good Point! I added two sentences and links for further info for now.
@@ -0,0 +1,1155 @@ | |||
{ |
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.
@@ -0,0 +1,1155 @@ | |||
{ |
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.
The comments here and in the next few cells could rather become markdown cells that guide through the process:
Minimize energy and output the minimized system
simulation.minimizeEnergy() with open(DATA / "topology.pdb", "w") as pdb_file: app.PDBFile.writeFile(simulation.topology, simulation.context.getState(getPositions=True, enforcePeriodicBox=True).getPositions(), file=pdb_file, keepIds=True)
Reply via ReviewNB
@@ -0,0 +1,1155 @@ | |||
{ |
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.
Is the number for the 'next' notebook already fixed? If yes, please include it here.
... refer to Txxx ... (maybe even with a link if available)
At the very and, we usually have a small quiz (three questions) can you think of any?
Reply via ReviewNB
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.
The number is T020.
We could link to the right folder in GitHub or do we wait for the website?
@Mika-Le the talktorial is GREAT!!! Very easy to follow and very good explanations, well done! |
Hi @AndreaVolkamer, thank you for the great feedback. I updated the notebook, it is ready to be reviewed again. |
@Mika-Le very well done (@schallerdavid), the notebook is great! I only added tiny textual adaptions. Unfortunately, I couldn't get it running on colab due to the mentioned As soon as the colab thing is solved, the notebook is ready to be merged! |
As per our chat, the issue is solved but we can't make sure it works reliably right now because Anaconda.org servers are a bit irresponsive now and installation fails. We have decided to reconvene in a week, since by then openmmforcefields might be on conda forge already (and does not suffer from server issues since it runs on separate network). |
|
@jaimergp, |
|
||
## Categories | ||
|
||
This talktorial is part of the following categories: [Collections overview](link) |
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 think this is outdated. Are we writing the categories in the notebook now?
@@ -0,0 +1,1250 @@ | |||
{ |
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.
Note that !command
does not stop the cell from being executed. "Dependencies successfully installed" will be showed regardless the success of the conda install
command. This, in addition to the cleared outputs, can be very misleading... We should check whether things were actually installed or not, maybe try importing openmm or something similar. Only if that succeeds we clear outputs and report success. Otherwise we want the logs to see what went wrong (e.g. the random network errors we saw).
Reply via ReviewNB
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.
Is it pythonic/ok to nest try except blocks? Something along the lines of
try: import condacolab [...] try: import rdkit clear_output() print("Installation successful") except: print("Dependencies not installed") except: print("Not on Colab")
If not, what would be the better way to check for successful installation? Any If-else-checks
possible for installed dependencies?
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'd use the else
clause of the try
block, which will only be reached if no exception was raised during try
try:
import condacolab
...
except:
print("Not on Colab")
else:
try:
import rdkit
clear_output()
print("Installation successful")
except:
print("Dependencies not installed!")
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.
Thanks, that is way better
@@ -0,0 +1,1250 @@ | |||
{ |
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.
@@ -0,0 +1,1250 @@ | |||
{ |
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.
RCSB offers a simple way to get PDBs via URLs: https://files.rcsb.org/download/3POZ.pdb
. Do we need the get_pdb_file
function? We could use requests
instead. Up to you!
Reply via ReviewNB
@@ -0,0 +1,1250 @@ | |||
{ |
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.
PDBFixer should be able to deal with Path objects but it does not now. We should open an issue for that.
Also, can we parameterize the pH value with an optional keyword? (..., pH=7.0):
Reply via ReviewNB
@@ -0,0 +1,1250 @@ | |||
{ |
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.
@@ -0,0 +1,1250 @@ | |||
{ |
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 thought OpenForceField could convert to OpenMM just fine. Am I missing something?
Reply via ReviewNB
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.
Since we have to convert the unit of the positions we extract topology and positions separately and create a modeller object from there. Maybe a modeller object could be created from the openff Molecule directly, but what about the unit conversion then? We could probably do it at another step of the process, if that's better. Any ideas?
@@ -0,0 +1,1250 @@ | |||
{ |
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.
Are you sure you need that +1 in complex_positions[len(protein.positions) + 1: total_atoms]
? Python indexing is [start, end)
, so I would assume that we can start the 2nd range where the 1st ended, not shifted +1. What's the extra atom? You can also drop 0
and total_atoms
. In other words:
complex_positions = unit.Quantity(np.zeros([total_atoms, 3]), unit=unit.nanometers) complex_positions[:len(protein.positions)] = protein.positions # add protein positions complex_positions[len(protein.positions):] = ligand.positions # add ligand positions
Reply via ReviewNB
@@ -0,0 +1,1250 @@ | |||
{ |
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.
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.
Didn't find anything suitable about tip3p, especially not in context with amber. Added a link to the wikipedia article about water models which includes a section for tip3p. Better recommendations are welcome
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 think this is the canonical citation: https://aip.scitation.org/doi/10.1063/1.445869
@@ -0,0 +1,1250 @@ | |||
{ |
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.
This # NBVAL_CHECK_OUTPUT
will almost certainly fail because the minimization is a stochastic process and we are drawing velocities from a distribution, so for sure we will not get the same exact numbers every time.
Remove this comment and maybe, in a new cell, check if everything ran correctly by checking that the simulation objects has reached the desired number of steps.
Reply via ReviewNB
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.
How about removing the #NBVAL_CHECK_OUTPUT from the simulation step cell and adding another cell afterwards with the following content:
simulation_time = round(simulation.context.getState().getTime().value_in_unit(unit.femtoseconds))
if simulation_time == steps * 2:
print("Simulation finished successfully!")
else:
print("Simulation failed!")
# NBVAL_CHECK_OUTPUT
This is awesome! I am sure OpenMM users will like this material a lot! Thanks for putting it together! I have added some minor comments that do not change the outcome or purpose of the talktorial, just some technical details we need to take into account. |
Ping me again when the comments are addressed so I resolve the conflicts and merge! |
Thanks @Mika-Le! 🚀 |
Details
Content review
here
.DataFrames
)Code review
a_variable_name
vsaVariableName
)black -l 99
)for i in range(len(list))
(see slides)# TODO: CI
import ...
lines are at the top (practice part) cell, ordered by standard library / 3rd party packages / our own (teachopencadd.*
)