Elastic network #54

jbarnoud · 2018-02-14T07:47:55Z

A first go at the rubber band elastic network. I have to figure out why I have discrepancies in my distance matrix compared to martinize1.

This is a step forward for #33. It requires #41 and #53 to be merged.

jbarnoud · 2018-02-14T10:03:33Z

The discrepancy in the distance matrix is likely due to #36. Because we do not have weights in the mapping, yet, the BB bead positions are wrong by a bit. Hence the distance matrix does not match martinize1.

jbarnoud · 2018-02-14T22:03:21Z

There are discrepancies in the mappings between martinize1 and backward. This, plus the fact that I was not using the masses to compute the centre of mass -_- was causing the errors in the elastic network since the beads were not located where they should have.

I also calculate the weights from the mapping as I thought it was the problem before realizing the mapping discrepancy. So this PR fixes #36.

Now, the elastic network is good, and the only differences are due to rounding. Note that it can lead to bonds being added or removed compare to martinize1 when the bond length is very close to the cutoff.

jbarnoud · 2018-02-16T12:33:33Z

Rebased against master after renaming.

pckroon · 2018-03-01T18:48:22Z

Also for this one, I'll review it once the other two are merged to prevent duplication of effort.

jbarnoud · 2018-03-12T19:14:14Z

Rebased against master.

* Add a `AttachMass` processor. * Compute mapping weights when reading the mapping files. The weights are attached to the bead node, alongside the underlying graphs. * `AverageBeads` uses the mapping-weighted COM instead of the unweighted COG

Histidine is still different than martinize1 because of an hydrogen. This should be fixed with protonation as PTMs.

Bonds are excluded for residue separated by 2 or less residues. The previous try only accounted for separations of 1.

During the mapping step, attributes from the original molecules are passed on to the mapped ones. This leads to the mapped molecules to have wrong attributes, such as atomistic masses on a coarse grained bead. This commit allows to keep only the attributes listed as argument. The attributes required to identify a residue are passed on in any case. These default attributes are chain, resid, and resname.

Some node attributes may be set to None when they are unknown, instead of being undefined. This means that `node.get(attr, 0)` can return `None` even if there is no known valid value. This breaks the PDB writer is some situations. This commit makes sure that the PDB writer uses its default value wjen an attribute is set to `None`.

jbarnoud · 2018-03-13T00:06:34Z

Rebased (again) against master. I'll be less responsive for 8 hours or so.

pckroon

Some minor things. In addition, there's docstrings missing, and I expect well find a dragon in do_mapping later on. Can't find it now though.

pckroon · 2018-03-13T01:31:35Z

vermouth/processors/apply_rubber_band.py

+
+
+def compute_decay(distance_matrix, shift, rate, power):
+ return np.exp(-rate * ((distance_matrix - shift) ** power))


I would call it just distance instead of distance_matrix, and let numpy deal with the broadcasting

pckroon · 2018-03-13T01:32:12Z

vermouth/processors/apply_rubber_band.py

+ decay_factor, decay_power, base_constant,
+ minimum_force):
+ constants = compute_decay(distance_matrix, lower_bound, decay_factor, decay_power)
+ np.fill_diagonal(constants, 0)


But do note that this should be distance_matrix here

I don't think so. The point of that line is to make sure that the force constant for the self pair is 0, which could not be the case depending of the parameters for the decay; but may especially not be the case because of rounding.

I meant at line 33 ;)

Oh! It was related to the comment above!

Jup :) You're too fast

pckroon · 2018-03-13T01:33:49Z

vermouth/processors/apply_rubber_band.py

+ np.fill_diagonal(constants, 0)
+ constants *= base_constant
+ constants[constants < minimum_force] = 0
+ constants[distance_matrix > upper_bound] = 0


constants[(constants < minimun_force) | (distance_matrix > upper_bound)] = 0 is also an option. Yours is more readable. Do you think it's worth profiling?

It would work indeed. I'll let optimization for latter, though. We need a general profiling pass anyway.

pckroon · 2018-03-13T01:41:41Z

vermouth/processors/apply_rubber_band.py

+ if separation <= 0:
+ raise ValueError('Separation has to be strictly positive.')
+ if separation == 1:
+ return nx.to_numpy_matrix(graph, nodelist=selection).astype(bool)


Could use a comment

pckroon · 2018-03-13T01:46:55Z

vermouth/processors/apply_rubber_band.py

+ subgraph = graph.subgraph(selection)
+ connectivity = np.zeros((len(subgraph), len(subgraph)), dtype=bool)
+ for (i, key_i), (j, key_j) in itertools.combinations(enumerate(subgraph.nodes), 2):
+ shortest_path = len(nx.shortest_path(subgraph, key_i, key_j))


I would've done it the other way around: make a complete distance matrix, and do numpy magic to find the node_keys between which the distance is what you want. There's not really a convenient method in networkx for this though. Maybe networkx.algorithms.shortest_paths.dense.floyd_warshall_numpy, but it would require profiling.

pckroon · 2018-03-13T01:54:34Z

vermouth/processors/apply_rubber_band.py

+ upper_bound, decay_factor, decay_power,
+ base_constant, minimum_force)
+ connectivity = build_connectivity_matrix(molecule, 2, selection=selection)
+ constants *= ~connectivity


This one needs a comment

pckroon · 2018-03-13T01:55:55Z

vermouth/processors/apply_rubber_band.py

+ from_key = selection[from_idx]
+ to_key = selection[to_idx]
+ force_constant = constants[from_idx, to_idx]
+ length = round(distance_matrix[from_idx, to_idx], 5)


Why not do distance_matrix = distance_matrix.round(5) before the loop?

pckroon · 2018-03-13T01:56:34Z

vermouth/processors/apply_rubber_band.py

+ to_key = selection[to_idx]
+ force_constant = constants[from_idx, to_idx]
+ length = round(distance_matrix[from_idx, to_idx], 5)
+ if force_constant > minimum_force:


Should this check not be part of compute_force_constants?

compute_force_constants compute the constants, what we do with them is not its responsibility. Also, filtering the force constants in compute_force_constants will only change the condition of the test, it will not get rid of it.

In make_force_constants there's already a check that does
constants[constants < minimum_force] = 0
So maybe just if force_constant?

Indeed. It is possible, but then the way it is here may be more explicit.

Hmmn. Or set it to np.nan in compute_force_constants, and check for that here. Now you're doing the same thing twice. And set a comment that it's set to nan in said function

pckroon · 2018-03-13T01:59:25Z

vermouth/processors/do_mapping.py

+ if from_idxs:
+ self.weights[to_idx][from_idxs[0]] = (
+ weights[(res_to, name_to)][(res_from, name_from)]
+ )


Why only from_idxs[0]? I'm lost in my own code :s

I am not completely sure either from the top of my head. But I'd rather have some tests in place before messing around with some refactoring.

jbarnoud

I guess I'm ought to do a massive run of docstringing soon.

jbarnoud · 2018-03-13T08:53:33Z

vermouth/processors/apply_rubber_band.py

+ decay_factor, decay_power, base_constant,
+ minimum_force):
+ constants = compute_decay(distance_matrix, lower_bound, decay_factor, decay_power)
+ np.fill_diagonal(constants, 0)


I don't think so. The point of that line is to make sure that the force constant for the self pair is 0, which could not be the case depending of the parameters for the decay; but may especially not be the case because of rounding.

jbarnoud · 2018-03-13T08:53:45Z

vermouth/processors/apply_rubber_band.py

+ np.fill_diagonal(constants, 0)
+ constants *= base_constant
+ constants[constants < minimum_force] = 0
+ constants[distance_matrix > upper_bound] = 0


It would work indeed. I'll let optimization for latter, though. We need a general profiling pass anyway.

jbarnoud · 2018-03-13T08:53:56Z

vermouth/processors/apply_rubber_band.py

+ to_key = selection[to_idx]
+ force_constant = constants[from_idx, to_idx]
+ length = round(distance_matrix[from_idx, to_idx], 5)
+ if force_constant > minimum_force:


compute_force_constants compute the constants, what we do with them is not its responsibility. Also, filtering the force constants in compute_force_constants will only change the condition of the test, it will not get rid of it.

jbarnoud · 2018-03-13T08:54:14Z

vermouth/processors/do_mapping.py

+ if from_idxs:
+ self.weights[to_idx][from_idxs[0]] = (
+ weights[(res_to, name_to)][(res_from, name_from)]
+ )


I am not completely sure either from the top of my head. But I'd rather have some tests in place before messing around with some refactoring.

pckroon · 2018-03-13T19:20:20Z

By the way, is the WIP label/tag still relevant?

jbarnoud · 2018-03-13T19:40:38Z

Hum... In a sense... everything is always a work in progress?

jbarnoud · 2018-03-20T10:59:04Z

What is still blocking in this PR?

pckroon · 2018-03-20T17:36:46Z

IIRC the missing docstrings.

jbarnoud · 2018-03-26T15:52:01Z

I added docstring!

pckroon

Excellent work! 2 minor comments and it's good to go. If you disagree, feel free to merge as is.

pckroon · 2018-03-26T17:34:48Z

vermouth/processors/apply_rubber_band.py



 def compute_decay(distance, shift, rate, power):
+ """
+ Compute the decay function of the force constant as function to the distance.


too vague imho

I will see with @Tsjerk if he can explain the parameters to me, or I will try to plot the function... This is an obscure feature that I never seen used.

It's already better if you describe the functional form. No need for shiny.

Well, little need.

pckroon · 2018-03-26T17:37:56Z

vermouth/processors/apply_rubber_band.py

+ The base force constant for the bonds in kJ/mol/nm^2. If 'decay_factor'
+ or 'decay_power' is set to 0, then it will be the used force constant.
+ minimum_force: float
+ Minimum force constat in kJ/mol/nm^2 under which bonds are not kept.


:math:? Or TeX formatting for the units? I'm pretty sure sphinx/rest can deal with it somehow.

Units look ugly in latex. It would be something like $$\frac{\text{kJ}}{\text{mol}\cdot\text{nm}^2}$$, it renders beautifully in latex (not in github it seems, but still), but it is super ugly in plain text.

I would probably write something liek this: $$\frac{kJ}{nm}^2$$; but I take your point. Is there no middle ground?

jbarnoud · 2018-03-26T17:53:05Z

I applied a compromise on the units.

jbarnoud added the WIP label Feb 14, 2018

jbarnoud added this to the Reproduce functionality from current Martinize milestone Feb 14, 2018

jbarnoud added this to To Do in Beta release via automation Feb 14, 2018

jbarnoud moved this from To Do to In progress in Beta release Feb 14, 2018

jbarnoud force-pushed the elastic branch from 6198fce to b2169ed Compare February 14, 2018 21:57

jbarnoud force-pushed the elastic branch from b2169ed to 607fc78 Compare February 16, 2018 12:33

jbarnoud mentioned this pull request Mar 7, 2018

Reworked parts of repair_graph to detect modifications #2

Closed

jbarnoud force-pushed the elastic branch from 607fc78 to fc756b3 Compare March 12, 2018 19:09

jbarnoud added 6 commits March 13, 2018 01:05

Adjust the mapping files to match martini1

b925124

Histidine is still different than martinize1 because of an hydrogen. This should be fixed with protonation as PTMs.

First go at the rubber band elastic network

3dd2d36

Fix connectivity in rubber band

9600bff

Bonds are excluded for residue separated by 2 or less residues. The previous try only accounted for separations of 1.

jbarnoud force-pushed the elastic branch from fc756b3 to d6dbcff Compare March 13, 2018 00:05

pckroon requested changes Mar 13, 2018

View reviewed changes

jbarnoud commented Mar 13, 2018

View reviewed changes

Address @pckroon's comments on #54

22858e4

jbarnoud removed the WIP label Mar 13, 2018

jbarnoud changed the title ~~[WIP] Elastic network~~ Elastic network Mar 13, 2018

Add documentation on the apply_rubber_band module

fd4510f

pckroon approved these changes Mar 26, 2018

View reviewed changes

Adress comments from @pckroon on docstrings

b286c84

pckroon approved these changes Mar 26, 2018

View reviewed changes

pckroon merged commit 0431abf into master Mar 26, 2018

Beta release automation moved this from In progress to Done Mar 26, 2018

jbarnoud deleted the elastic branch March 28, 2018 09:07

pckroon mentioned this pull request Nov 2, 2023

-el fix #559

Merged



		def compute_decay(distance_matrix, shift, rate, power):
		return np.exp(-rate * ((distance_matrix - shift) ** power))

Elastic network #54

Elastic network #54

Conversation

jbarnoud commented Feb 14, 2018

jbarnoud commented Feb 14, 2018

jbarnoud commented Feb 14, 2018

jbarnoud commented Feb 16, 2018

pckroon commented Mar 1, 2018

jbarnoud commented Mar 12, 2018

jbarnoud commented Mar 13, 2018

pckroon left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pckroon Mar 13, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jbarnoud left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pckroon commented Mar 13, 2018

jbarnoud commented Mar 13, 2018

jbarnoud commented Mar 20, 2018

pckroon commented Mar 20, 2018

jbarnoud commented Mar 26, 2018

pckroon left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jbarnoud commented Mar 26, 2018

pckroon Mar 13, 2018 •

edited

Loading