Skip to content

Commit

Permalink
Merge branch 'clean_flat' into 'master'
Browse files Browse the repository at this point in the history 
Clean flat

See merge request santuz/heligeom_webserver!10
  • Loading branch information
Hubert Santuz committed Dec 12, 2024
2 parents 99fb4bc + 9e8c0bc commit 500bded
Show file tree
Hide file tree
Showing 9 changed files with 584 additions and 491 deletions.
2 changes: 1 addition & 1 deletion .gitmodules
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
[submodule "ptools-python"]
path = ptools-python
url = https://github.com/ptools/ptools.git
branch = develop
branch = main
851 changes: 446 additions & 405 deletions Pipfile.lock
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion heligeom/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ def add_examples(db):
db.session.commit()


def create_app(config_object=None, clear_database=True):
def create_app(config_object=None, clear_database=False):
# create and configure the app
app = Flask(__name__, instance_relative_config=True)
app.config.from_object(config_object)
Expand Down
4 changes: 2 additions & 2 deletions heligeom/templates/macros/forms.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,14 +42,14 @@
{{ chain }}
<label>Chain ID</label>
{% if chain.errors %}
<span class="helper-text error-input">{{ chain.errors[0] }}</span>
<span class="helper-text error-input">{{ chain.errors }}</span>
{% endif %}
</div>
<div class="col s12 m6 input-field">
{{ res_range }}
<label>Residue Range</label>
{% if res_range.errors %}
<span class="helper-text error-input">{{ res_range.errors[0] }}</span>
<span class="helper-text error-input">{{ res_range.errors }}</span>
{% endif %}
</div>
</div>
Expand Down
16 changes: 16 additions & 0 deletions heligeom/templates/results_2_oligomers.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,14 @@ <h4 class="header center">2nd Oligomeric form</h4>
});
</script>
</div>
{% if construct_data.flatten %}
<div class="row">
<h5 class="left-align">Results of the flatten algorithm</h5>
<div class="col s12 m6">
{{ results_macro.display_ring_params(construct_data) }}
</div>
</div>
{% endif %}
</div> <!-- end card content-->
</div> <!-- end card -->
</div> <!-- end col -->
Expand All @@ -191,6 +199,14 @@ <h4 class="header center">2nd Oligomeric form</h4>
});
</script>
</div>
{% if construct_data_bis.flatten %}
<div class="row">
<h5 class="left-align">Results of the flatten algorithm</h5>
<div class="col s12 m6">
{{ results_macro.display_ring_params(construct_data_bis) }}
</div>
</div>
{% endif %}
</div> <!-- end card content-->
</div> <!-- end card -->
</div> <!-- end col -->
Expand Down
2 changes: 1 addition & 1 deletion heligeom/tool.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,7 +315,7 @@ def heli_ring(self, fileout, ncopies, z_align=True):
monomers_per_turn = round(360.0 / abs(math.degrees(self.hp.angle)))

# Call adjust with Ptarget = 0 to create a ring
newhp, newARb, bestener, rms, fnat = utils.adjust(arec, self.hp, eref, monomers_per_turn, 0)
newhp, rms, fnat = utils.adjust(arec, self.hp, eref, monomers_per_turn, 0)

self.oligomer = heli_construct(self.monomer1.rb, newhp, N=ncopies, Z=z_align)
io.write_pdb(self.oligomer, fileout) # type: ignore
Expand Down
175 changes: 97 additions & 78 deletions heligeom/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,11 +87,54 @@ def check_pair_monomers(chain1, res_range1, chain2, res_range2):


def ener1(rec, lig, cutoff):
"""Compute the energy between 2 AttractRigidBody from an AttractForceField.
Parameters
----------
rec : AttractRigidBody
1st element
lig : AttractRigidBody
2nd element
cutoff : float
cut off for non bonded terms.
Returns
-------
float
the computed energy
"""
ff = forcefield.AttractForceField1(rec, lig, cutoff)
return ff.energy()


def adjust(ARb, hpori, enref, Ntarget, Ptarget):
"""Adjust the Screw object to match the Ntarget & Ptarget.
The algorithm use 2 MonteCarlo cycle to find the best solution.
Parameters
----------
ARb : AttractRigidBody
_description_
hpori : Screw
the original Screw to modify
enref : float
Reference energy
Ntarget : float
The number of monomers per turn to achieve
Ptarget : float
The pitch to achieve
Returns
-------
Screw
the new screw object
Float
the RMSD between the original Rb and the new one
Float
The fNaT between the original Rb and the new one
"""
rec = ARb.copy()
p = measure.center(rec)

Expand All @@ -107,7 +150,6 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
transform.translate(ligor2, -1 * axe * hpori.normtranslation)

# computes distance from the axis, rayon, radial axis vectn
# rayori = math.sqrt((pt0.x-p.x)**2 + (pt0.y-p.y)**2 + (pt0.z-p.z)**2) # Hub: unused
vect = pt0 - p
vectn = vect - np.dot(vect, axe) * axe
rayon = np.linalg.norm(vectn)
Expand All @@ -125,9 +167,7 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):

# computes new distance from the axis
raynew = rayon * math.sin(hpori.angle / 2) / math.sin(angnew / 2)
pt0new = p + raynew * vectn # --> axis is laterally translated???
print("angini,transini", hpori.angle, hpori.normtranslation)
print("angnew,transnew,rayon", angnew, transnew, raynew)
pt0new = p + raynew * vectn

##
lig = rec.copy()
Expand All @@ -136,8 +176,6 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
transform.translate(lig, axe * transnew)
eninit = ener1(rec, lig, 7)

# debug
# print "\nAfter modifying ang and trans :\n"
print(
"#INIT ",
Ptarget,
Expand All @@ -148,53 +186,51 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
f" - RMSD ligor ligTarget: {min(rmsd(ligor1, lig), rmsd(ligor2, lig)):7.2f}",
)
print(
f"# Fnat rec ligTarget: {max(measure.fnat(rec, ligor1, rec, lig, 7), measure.fnat(rec, ligor2, rec, lig, 7)):7.2f}"
"#Fnat rec ligTarget: ",
f"{max(measure.fnat(rec, ligor1, rec, lig, 7), measure.fnat(rec, ligor2, rec, lig, 7)):7.2f}",
)
# debug

nbiter = 300
nbiter2 = 600
# MC iterations
iter_1stMC = 300
iter_2ndMC = 600
best_iter = 0
nbaccept = 0
nbtry = 0

# variables for random generation
uang = 5.0
udis = 3.0
gm = 0.0
gsda = 2.0
gsdd = 0.5
print()
print("Parameters: ", nbiter, uang, udis, "\t", nbiter2, gm, gsda, gsdd)
print()
dax, day, daz = 0.0, 0.0, 0.0

# nbiter=500
# nbaccept=0.
raytmp = rayon
enold = eninit
# Init variables for MC

pt0tmp = pt0.copy()
dax, day, daz = 0.0, 0.0, 0.0
# daxbst, daybst, dazbst = 0., 0., 0. # Hub unused
bestene = 9999
enold = eninit

raybst = raynew
raytmp = rayon

pt0bst = np.array([0.0, 0.0, 0.0])
pt0tmp = pt0.copy()

recbst = rec.copy()
rectmp = rec.copy()
recnew = rec.copy()
fnatbst = measure.fnat(rec, ligor1, rec, lig, 7) # Hub : 7 as cutoff?

# MC search - constant (N,P) ==> (angnew,transnew)
fnatbst = measure.fnat(rec, ligor1, rec, lig, 7)

ibst = 0
nbaccept = 0
nbtry = 0
for i in range(nbiter):
# MC search - constant (N,P) ==> (angnew,transnew)
for i in range(iter_1stMC):
if i % 100 == 0:
print(f"... {i}")
rectmp = recnew.copy()
if random.random() > 0.5: # rotation rec; lig will follow
dax = math.radians(random.uniform(-1 * uang, uang)) # gaussian rather than uniform???
dax = math.radians(random.uniform(-1 * uang, uang))
day = math.radians(random.uniform(-1 * uang, uang))
daz = math.radians(random.uniform(-1 * uang, uang))
transform.ab_rotate(
rectmp, p, p + vectn, dax
) # p = pinit... mais c'est l'axe qui bouge...
transform.ab_rotate(rectmp, p, p + vectn, dax, degrees=False)
transform.ab_rotate(rectmp, p, p + wectn, day, degrees=False)
transform.ab_rotate(rectmp, p, p + axe, daz, degrees=False)
else: # translation rec; lig will follow
Expand All @@ -206,12 +242,10 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
transform.ab_rotate(ligtmp, pt0tmp, pt0tmp + axe, angnew, degrees=False)
transform.translate(ligtmp, axe * transnew)

# ICI TEST FNAT: si FNAT < 0.25 on passe le tour
# attention, rec a bouge, il faut utiliser fnat2
fn = measure.fnat(rec, rectmp, ligor1, ligtmp, 7) # Hub : cutoff = 7
if fn < 0.3:
print("fnat too low ", i, dax, day, daz, raynew - raytmp, fn)
fnat = measure.fnat(rec, rectmp, ligor1, ligtmp, 7)
if fnat < 0.3:
continue # goto next iteration

nbtry += 1
enew = ener1(rectmp, ligtmp, 7)

Expand All @@ -224,8 +258,8 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
pt0bst = pt0tmp
recbst = rectmp.copy()
raybst = raytmp
fnatbst = fn
ibst = i
fnatbst = fnat
best_iter = i
if random.random() <= math.exp(-delta): # always true if delta == 0
enold = enew
recnew = rectmp
Expand All @@ -234,44 +268,38 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):

# construct final model
# ligand
bestmp = recbst.copy()
transform.ab_rotate(bestmp, pt0bst, pt0bst + axe, angnew, degrees=False)
transform.translate(bestmp, axe * transnew)
# superpose recbst on rec and applymatrix on final and best [CP 12.06.24: superpose the axes instead?]
ligbst = bestmp.copy()
matfi = superpose.fit_matrix(rec, recbst)
transform.move(ligbst, matfi)
ligbst = recbst.copy()
transform.ab_rotate(ligbst, pt0bst, pt0bst + axe, angnew, degrees=False)
transform.translate(ligbst, axe * transnew)

# hpnew = mat_trans_to_screw(superpose(rec,ligbst).matrix) CORR 4.11.22
hpbst = heli_analyze(recbst, ligbst)

# debug
print()
print("#ADJ Nbest raynew raybst-rayinit bestnrj rmsd(ligori,ligbst) fnat")
min_rmsd = min(rmsd(ligor1, ligbst), rmsd(ligor2, ligbst))
print(
f"#ADJ2 ${ibst:3d} ${raybst:10.2f} ${raybst-rayon:15.2f} {bestene:12.2f} {min_rmsd:15.2f} ${fnatbst:12.2f}"
"#ADJ2 ",
f"{best_iter:3d} {raybst:10.2f} {raybst-rayon:15.2f} {bestene:12.2f} ",
f"{min_rmsd:15.2f} {fnatbst:12.2f}",
)
print(f"acc: {nbaccept / nbtry} {nbtry}")
# debug

print("acc: ", nbaccept / nbtry, nbtry)

# new MC run around the best energy
recnew = recbst
recnew = recbst.copy()
raynew = raybst
nbtry2 = 0.0
nbaccept2 = 0.0
for i in range(nbiter2):
for i in range(iter_2ndMC):
if i % 100 == 0:
print(f"... {i}")
rectmp = recnew.copy()
if random.random() > 0.5: # rotation rec; lig will follow
dax = math.radians(random.gauss(0.0, gsda)) # gaussian rather than uniform???
dax = math.radians(random.gauss(0.0, gsda))
day = math.radians(random.gauss(0.0, gsda))
daz = math.radians(random.gauss(0.0, gsda))
transform.ab_rotate(
rectmp, p, p + vectn, dax
) # p = pinit... mais c'est l'axe qui bouge...
transform.ab_rotate(rectmp, p, p + vectn, dax, degrees=False)
transform.ab_rotate(rectmp, p, p + wectn, day, degrees=False)
transform.ab_rotate(rectmp, p, p + axe, daz, degrees=False)
else: # translation rec; lig will follow
Expand All @@ -283,19 +311,12 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
transform.ab_rotate(ligtmp, pt0tmp, pt0tmp + axe, angnew, degrees=False)
transform.translate(ligtmp, axe * transnew)

# ICI TEST FNAT: si FNAT < 0.25 on passe le tour
# attention, rec a bouge, il faut utiliser fnat2
fn = measure.fnat(rec, rectmp, ligor1, ligtmp, 7) # Hub : cutoff = 7
if fn < 0.3:
print("fnat too low ", i, dax, day, daz, raynew - raytmp, fn)
fnat = measure.fnat(rec, rectmp, ligor1, ligtmp, 7)
if fnat < 0.3:
continue # goto next iteration
nbtry2 += 1
enew = ener1(rectmp, ligtmp, 7)

# DEGUG
if enew < -170:
print("i, enew: ", i, enew)

# MC test
delta = enew - enold
if delta < 0:
Expand All @@ -305,8 +326,8 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):
pt0bst = pt0tmp
recbst = rectmp.copy()
raybst = raytmp
fnatbst = fn
ibst = i
fnatbst = fnat
best_iter = i
if random.random() <= math.exp(-delta): # always true if delta == 0
enold = enew
recnew = rectmp
Expand All @@ -315,32 +336,30 @@ def adjust(ARb, hpori, enref, Ntarget, Ptarget):

# construct final model
# ligand
bestmp = recbst.copy()
transform.ab_rotate(bestmp, pt0bst, pt0bst + axe, angnew, degrees=False)
transform.translate(bestmp, axe * transnew)
ligbst = recbst.copy()
transform.ab_rotate(ligbst, pt0bst, pt0bst + axe, angnew, degrees=False)
transform.translate(ligbst, axe * transnew)

# superpose recbst on rec and applymatrix on final and best [CP 12.06.24: superpose the axes instead?]
ligbst = bestmp.copy()
matfi = superpose.fit_matrix(rec, recbst)
matfi = superpose.fit_matrix(recbst, rec)
transform.move(ligbst, matfi)

# hpnew = mat_trans_to_screw(superpose(rec,ligbst).matrix) CORR 4.11.22
hpbst = heli_analyze(rec, ligbst)

# debug
print("#ADJ2 Nbest raynew raybst-rayinit bestnrj rmsd(ligori,ligbst) fnat")
min_rmsd = min(rmsd(ligor1, ligbst), rmsd(ligor2, ligbst))
print(
f"#ADJ2 ${ibst:3d} ${raybst:10.2f} ${raybst-rayon:15.2f} {bestene:12.2f} {min_rmsd:15.2f} ${fnatbst:12.2f}"
"#ADJ2 ",
f"{best_iter:3d} {raybst:10.2f} {raybst-rayon:15.2f} {bestene:12.2f} ",
f"{min_rmsd:15.2f} {fnatbst:12.2f}",
)
print(f"acc2: {nbaccept2 / nbtry2} {nbtry2}")
# debug

print("acc2: ", (nbaccept2 * 1.0) / (nbtry2 * 1.0), nbtry2)

monomers_per_turn = round(360.0 / abs(math.degrees(hpbst.angle)))
pitch = abs(hpbst.normtranslation * (360.0 / (abs(math.degrees(hpbst.angle)))))

print(f"hpbst : {hpbst.normtranslation}, {hpbst.angle}")
print(f"hpbst norm & angle : {hpbst.normtranslation}, {hpbst.angle}")
print(f"hpbst : mono, pitch: {monomers_per_turn}, {pitch} ")

return hpbst, ligbst, bestene, min(rmsd(ligor1, ligbst), rmsd(ligor2, ligbst)), fnatbst
return hpbst, min_rmsd, fnatbst
Loading

0 comments on commit 500bded

Please sign in to comment.