From 5dd9efe70c2b37b7797a1c09c4ab571195e54c62 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sat, 11 Nov 2023 22:38:48 -0700
Subject: [PATCH 01/17] Skeleton code.

---
 src/classes/constants.h |   5 ++
 src/classes/protein.cpp | 178 +++++-----------------------------------
 src/classes/protein.h   |   2 +-
 3 files changed, 28 insertions(+), 157 deletions(-)

diff --git a/src/classes/constants.h b/src/classes/constants.h
index a679c7a8b..622a8d2a6 100755
--- a/src/classes/constants.h
+++ b/src/classes/constants.h
@@ -318,6 +318,11 @@
 #define homology_long_axis_rotations 0
 #define homology_region_optimization 0
 
+#define reconnect_length_limit 5
+#define reconnect_generations 8
+#define reconnect_keepbest 4
+#define reconnect_angle_importance 0.1
+
 //////////////////////////////////////////////////////////////////////////////////////////////
 
 // Debugging stuff.
diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index 0e030f167..2774eba88 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1815,168 +1815,34 @@ void Protein::conform_backbone(int startres, int endres,
     mass_undoable = wmu;
 }
 
-#define DBG_BCKCONN 0
-#define _INCREMENTAL_BKCONN 1
-void Protein::backconnect(int startres, int endres)
+float Protein::reconnect(int startres, int endres)
 {
-    save_undo_state();
-    bool wmu = mass_undoable;
-    mass_undoable = true;
-
-    int i;
-    int inc = sgn(endres-startres);
-
-    #if _INCREMENTAL_BKCONN
-    int iter;
-    for (iter=24; iter>=0; iter--)
+    int length = abs(endres - startres) + 1;
+    if (length > reconnect_length_limit)
     {
-    #endif
-
-        // Glue the last residue onto the target.
-        // Then adjust its inner bonds so the other end points as closely to the previous residue as possible.
-        // Then do the same for the previous residue, and the one before, etc,
-        // all the way back to the starting residue.
-        // Give a warning if the starting residue has an anomaly > 0.1A.
-        AminoAcid *next, *curr, *prev;
-        int pointer = endres;
-        float movfactor = 1, decrement, anomaly = 0;
-
-        decrement = 1.0 / fabs(endres - startres);		// if exterior is the opposite of interior, what's the opposite of increment?
-
-        next = get_residue(endres+inc);
-        curr = get_residue(pointer);
-        prev = get_residue(pointer-inc);
-
-        #if DBG_BCKCONN
-        cout << "backconnect( " << startres << ", " << endres << ")" << endl;
-        #endif
-        while (next && curr)
-        {
-            #if DBG_BCKCONN
-            cout << pointer << ":";
-            #endif
-
-            Point pts[4];
-            (inc > 0) ? next->predict_previous_COCA(pts) : next->predict_next_NHCA(pts);
-            Point ptsc[4];
-
-            #if _INCREMENTAL_BKCONN
-            ptsc[0] = (inc > 0) ? curr->get_atom_location("C") : curr->get_atom_location("N");
-            ptsc[1] = (inc > 0) ? curr->get_atom_location("O") : curr->HN_or_substitute_location();
-            ptsc[2] = curr->get_atom_location("CA");
-
-            Point _4avg[3];
-            for (i=0; i<3; i++)
-            {
-                _4avg[0] = pts[i];
-                _4avg[1] = ptsc[i];
-                _4avg[0].weight = movfactor;
-                _4avg[1].weight = 1.0 - movfactor;
-                pts[i] = average_of_points(_4avg, 2);
-            }
-            #endif
-
-            curr->attach_to_prediction(pts, inc > 0);
-            // break;
-
-            #if DBG_BCKCONN
-            cout << "g";
-            #endif
-
-            if (prev)
-            {
-                MovabilityType fmov = curr->movability;
-                curr->movability = MOV_ALL;
-
-                #if DBG_BCKCONN
-                cout << "a";
-                #endif
-
-                pts[4];
-                (inc < 0) ? prev->predict_previous_COCA(pts) : prev->predict_next_NHCA(pts);
-                Point target_heavy = pts[0];
-                Point target_pole = pts[1];
-
-                #if DBG_BCKCONN
-                cout << "b";
-                #endif
-
-                float theta, step = fiftyseventh*1.0, r, btheta = 0, bestr;
-                for (theta=0; theta < M_PI*2; theta += step)
-                {
-                    r = target_heavy.get_3d_distance( (inc > 0) ? curr->get_atom_location("N") : curr->get_atom_location("C") );
-                    // r -= target_pole.get_3d_distance( (inc > 0) ? curr->HN_or_substitute_location() : curr->get_atom_location("O") );
-                    if (inc > 0) r -= target_pole.get_3d_distance(curr->HN_or_substitute_location());
-
-                    if (!theta || (r < bestr))
-                    {
-                        bestr = r;
-                        btheta = theta;
-                    }
-
-                    curr->rotate_backbone( (inc > 0) ? CA_desc : N_asc, step );
-                }
-                curr->rotate_backbone( (inc > 0) ? CA_desc : N_asc, btheta );
-
-                #if DBG_BCKCONN
-                cout << "c";
-                #endif
-
-                btheta=0;
-                for (theta=0; theta < M_PI*2; theta += step)
-                {
-                    r = target_heavy.get_3d_distance( (inc > 0) ? curr->get_atom_location("N") : curr->get_atom_location("C") );
-                    // r -= target_pole.get_3d_distance( (inc > 0) ? curr->HN_or_substitute_location() : curr->get_atom_location("O") );
-                    if (inc < 0) r -= target_pole.get_3d_distance(curr->get_atom_location("O"));
-
-                    if (!theta || (r < bestr))
-                    {
-                        bestr = r;
-                        btheta = theta;
-                    }
-
-                    curr->rotate_backbone( (inc > 0) ? C_desc : CA_asc, step );
-                }
-                curr->rotate_backbone( (inc > 0) ? C_desc : CA_asc, btheta );
-                anomaly = bestr;
-
-                #if DBG_BCKCONN
-                cout << "d";
-                #endif
-
-                curr->movability = fmov;
-            }
-
-            if (pointer == startres)
-            {
-                #if DBG_BCKCONN
-                cout << ". ";
-                #endif
-                break;
-            }
-
-            pointer -= inc;
-            next = curr;
-            curr = prev;
-            prev = get_residue(pointer-inc);
-            movfactor -= decrement;
-
-            #if DBG_BCKCONN
-            cout << ", ";
-            #endif
-        };
+        cout << "Attempt to reconnect strand in excess of length limit." << endl;
+        throw 0xb5712a9d;
+    }
+    int connres = (startres > endres) ? (endres - 1) : (endres + 1);
+    float angle = triangular;
+    int gen;
+    int max_candidates = pow(3, length) * reconnect_keepbest;
+    float candidates[reconnect_keepbest][max_candidates+1][length+1];
+    float best_candidates[reconnect_keepbest][length+1];
+    float best_score[reconnect_keepbest];
 
-        #if DBG_BCKCONN
-        cout << endl;
-        #endif
+    for (gen=0; gen<reconnect_generations; gen++)
+    {
+        // TODO: try every possible combination of starting phi,psi angles +/- the angle variable.
+        // Keep the best reconnect_keepbest from each generation and retry.
+        // Score each candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle.
 
-        /*if (anomaly > 0.1) cout << "Warning! conform_backbone( " << startres << ", " << endres << " ) anomaly out of range." << endl
-        						<< "# " << (startres-inc) << " anomaly: " << anomaly << endl;*/
-        #if _INCREMENTAL_BKCONN
+        angle /= 2;
     }
-        #endif
 
-    mass_undoable = wmu;
+    // TODO: Use the best candidate and trial-and-error it onto connres.
+
+    // TODO: Return the anomaly.
 }
 
 void Protein::make_helix(int startres, int endres, float phi, float psi)
diff --git a/src/classes/protein.h b/src/classes/protein.h
index 67bf84de8..5f529b335 100755
--- a/src/classes/protein.h
+++ b/src/classes/protein.h
@@ -155,7 +155,7 @@ class Protein
                           int iters, bool backbone_atoms_only
                          );
 
-    void backconnect(int startres, int endres);
+    float reconnect(int startres, int endres);
     void find_residue_initial_bindings();
     void undo();
 

From 280957e10cce7d86d928201e1a0ba8508bc8d604 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sat, 11 Nov 2023 22:53:56 -0700
Subject: [PATCH 02/17] Loops.

---
 src/classes/aminoacid.cpp |  6 ++++--
 src/classes/aminoacid.h   |  4 ++--
 src/classes/protein.cpp   | 36 ++++++++++++++++++++++++++++++++----
 3 files changed, 38 insertions(+), 8 deletions(-)

diff --git a/src/classes/aminoacid.cpp b/src/classes/aminoacid.cpp
index 6bb681f2f..1dc56fa1b 100755
--- a/src/classes/aminoacid.cpp
+++ b/src/classes/aminoacid.cpp
@@ -2171,23 +2171,25 @@ float AminoAcid::get_omega()
     return find_angle_along_vector(CA0->get_location(), CA1->get_location(), N->get_location(), axis);
 }
 
-bond_rotation_fail_reason AminoAcid::rotate_phi(float a)
+bond_rotation_fail_reason AminoAcid::rotate_phi(float a, bool absolute)
 {
     if (aadef->proline_like) return bf_disallowed_rotation;
     Atom* N  = get_atom("N");
     Atom* CA = get_atom("CA");
     Bond* b = CA->get_bond_between(N);
     if (!b) return bf_bond_not_found;
+    if (absolute) a -= get_phi();
     b->rotate(a, true, true);
     return b->last_fail;
 }
 
-bond_rotation_fail_reason AminoAcid::rotate_psi(float a)
+bond_rotation_fail_reason AminoAcid::rotate_psi(float a, bool absolute)
 {
     Atom* CA = get_atom("CA");
     Atom* C  = get_atom("C");
     Bond* b = CA->get_bond_between(C);
     if (!b) return bf_bond_not_found;
+    if (absolute) a -= get_psi();
     b->rotate(a, true, true);
     return b->last_fail;
 }
diff --git a/src/classes/aminoacid.h b/src/classes/aminoacid.h
index c97575658..aaea71942 100755
--- a/src/classes/aminoacid.h
+++ b/src/classes/aminoacid.h
@@ -148,8 +148,8 @@ class AminoAcid : public Molecule
     float get_phi();
     float get_psi();
     float get_omega();
-    bond_rotation_fail_reason rotate_phi(float rel_angle);
-    bond_rotation_fail_reason rotate_psi(float rel_angle);
+    bond_rotation_fail_reason rotate_phi(float rel_angle, bool absolute = false);
+    bond_rotation_fail_reason rotate_psi(float rel_angle, bool absolute = false);
 
     // Intermol functions.
     float get_intermol_binding(AminoAcid* neighbor, bool backbone_atoms_only = false);
diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index 2774eba88..9170593ed 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1823,19 +1823,47 @@ float Protein::reconnect(int startres, int endres)
         cout << "Attempt to reconnect strand in excess of length limit." << endl;
         throw 0xb5712a9d;
     }
-    int connres = (startres > endres) ? (endres - 1) : (endres + 1);
+    int dir = sgn(endres - startres);
+    int connres = (dir<0) ? (endres - 1) : (endres + 1);
     float angle = triangular;
     int gen;
-    int max_candidates = pow(3, length) * reconnect_keepbest;
+    int max_candidates = pow(3, length);
     float candidates[reconnect_keepbest][max_candidates+1][length+1];
     float best_candidates[reconnect_keepbest][length+1];
     float best_score[reconnect_keepbest];
 
+    int ikept, icand, ires, jres;
+    int working_candidate[length+1];
+
     for (gen=0; gen<reconnect_generations; gen++)
     {
         // TODO: try every possible combination of starting phi,psi angles +/- the angle variable.
-        // Keep the best reconnect_keepbest from each generation and retry.
-        // Score each candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle.
+        for (ikept=0; ikept < (gen ? reconnect_keepbest : 1); ikept++)
+        {
+            for (ires = 0; ires < length; ires++) working_candidate[ires] = -1;
+            for (icand = 0; icand < max_candidates; icand++)
+            {
+                for (ires = 0; ires < length; ires++)
+                {
+                    working_candidate[ires]++;
+                    if (working_candidate[ires] > 1) working_candidate[ires] = -1;
+                    else break;
+                }
+
+                for (ires = 0; ires < length; ires++)
+                {
+                    int resno = startres + dir*(ires>>1);
+                    candidates[ikept][icand][ires] = 
+                        ((ires & 0x1) ? get_residue(resno)->get_psi() : get_residue(resno)->get_phi())
+                        + angle * working_candidate[ires];
+
+                    // TODO: Rotate and score the backbone.
+                    // Score each candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle.
+                }
+
+                // TODO: Keep the best reconnect_keepbest from each generation and retry.
+            }
+        }
 
         angle /= 2;
     }

From a98f62b1a16228639426e7aaa9396260855264d0 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sat, 11 Nov 2023 23:04:04 -0700
Subject: [PATCH 03/17] Test/debug code.

---
 src/classes/protein.cpp | 5 ++++-
 src/interpreter.cpp     | 8 +++++---
 test/reconnect.pepd     | 8 ++++++++
 3 files changed, 17 insertions(+), 4 deletions(-)
 create mode 100644 test/reconnect.pepd

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index 9170593ed..4e68c3408 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1854,12 +1854,15 @@ float Protein::reconnect(int startres, int endres)
                 {
                     int resno = startres + dir*(ires>>1);
                     candidates[ikept][icand][ires] = 
-                        ((ires & 0x1) ? get_residue(resno)->get_psi() : get_residue(resno)->get_phi())
+                        ((ires & 0x1) ? get_residue(resno)->get_psi() : get_residue(resno)->get_phi())      // TODO: Persist previous best candidates.
                         + angle * working_candidate[ires];
 
+                    cout << candidates[ikept][icand][ires]*fiftyseven << " ";
+
                     // TODO: Rotate and score the backbone.
                     // Score each candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle.
                 }
+                cout << endl;
 
                 // TODO: Keep the best reconnect_keepbest from each generation and retry.
             }
diff --git a/src/interpreter.cpp b/src/interpreter.cpp
index 1864b3d04..6cdc012fd 100644
--- a/src/interpreter.cpp
+++ b/src/interpreter.cpp
@@ -1073,10 +1073,10 @@ int main(int argc, char** argv)
                 sr = interpret_single_int(words[l++]);
                 if (!words[l]) raise_error("Insufficient parameters given for CONNECT.");
                 ct = interpret_single_int(words[l++]);
-                if (words[l]) iters = interpret_single_int(words[l++]);
+                // if (words[l]) iters = interpret_single_int(words[l++]);
                 if (words[l]) raise_error("Too many parameters given for CONNECT.");
                 er = ct - sgn(ct-sr);
-
+/*
                 Atom *a1, *a2, *a3;
                 AminoAcid *cta = working->get_residue(ct), *era = working->get_residue(er);
 
@@ -1099,7 +1099,9 @@ int main(int argc, char** argv)
                 a3 = era->get_atom("CA");
 
                 working->conform_backbone(sr, er, a1, pt3[0], a2, pt3[1], iters);
-                working->backconnect(sr, er);
+                working->backconnect(sr, er);*/
+
+                working->reconnect(sr, er);
 
             _no_connect:
                 goto _pc_continue;
diff --git a/test/reconnect.pepd b/test/reconnect.pepd
new file mode 100644
index 000000000..d99fe6a7c
--- /dev/null
+++ b/test/reconnect.pepd
@@ -0,0 +1,8 @@
+
+LOAD "pdbs/OR51/OR51E2.upright.pdb"
+
+LET @axis = @5.50 - @6.48
+ECHO @axis
+
+ROTATE @axis -15 %6.48 %6.28 %6.59
+CONNECT %6.63 %6.59

From 83dac273d90101a4f76ab492ecaebcd0e7406753 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sun, 12 Nov 2023 07:47:16 -0700
Subject: [PATCH 04/17] Scoring.

---
 src/classes/protein.cpp | 65 +++++++++++++++++++++++++++++++++++++----
 1 file changed, 60 insertions(+), 5 deletions(-)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index 4e68c3408..5c3b8f9f5 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1823,18 +1823,54 @@ float Protein::reconnect(int startres, int endres)
         cout << "Attempt to reconnect strand in excess of length limit." << endl;
         throw 0xb5712a9d;
     }
-    int dir = sgn(endres - startres);
+
+    int dir = sgn(endres - startres); if (!dir) dir = 1;
+
     int connres = (dir<0) ? (endres - 1) : (endres + 1);
+    AminoAcid* connaa = get_residue(connres);
+    if (!connaa)
+    {
+        cout << "Attempt to reconnect strand to nonexistent target." << endl;
+        throw 0xb5712a9d;
+    }
+
+    AminoAcid* endaa = get_residue(endres);
+    if (!endaa)
+    {
+        cout << "Attempt to reconnect discontinuous strand." << endl;
+        throw 0xb5712a9d;
+    }
+
     float angle = triangular;
     int gen;
     int max_candidates = pow(3, length);
     float candidates[reconnect_keepbest][max_candidates+1][length+1];
     float best_candidates[reconnect_keepbest][length+1];
     float best_score[reconnect_keepbest];
+    float anomaly;
 
     int ikept, icand, ires, jres;
     int working_candidate[length+1];
 
+    Point should_be_at[4];
+    Atom* is_at[4];
+    if (dir > 0)
+    {
+        connaa->predict_previous_COCA(should_be_at);
+        is_at[0] = endaa->get_atom("C");
+        is_at[1] = endaa->get_atom("O");
+        is_at[2] = endaa->get_atom("CA");
+        is_at[3] = nullptr;
+    }
+    else
+    {
+        connaa->predict_next_NHCA(should_be_at);
+        is_at[0] = endaa->get_atom("N");
+        is_at[1] = endaa->HN_or_substitute();
+        is_at[2] = endaa->get_atom("CA");
+        is_at[3] = nullptr;
+    }
+
     for (gen=0; gen<reconnect_generations; gen++)
     {
         // TODO: try every possible combination of starting phi,psi angles +/- the angle variable.
@@ -1853,15 +1889,34 @@ float Protein::reconnect(int startres, int endres)
                 for (ires = 0; ires < length; ires++)
                 {
                     int resno = startres + dir*(ires>>1);
+                    AminoAcid* aa = get_residue(resno);
+                    if (!aa)
+                    {
+                        cout << "Attempt to reconnect discontinuous strand." << endl;
+                        throw 0xb5712a9d;
+                    }
                     candidates[ikept][icand][ires] = 
-                        ((ires & 0x1) ? get_residue(resno)->get_psi() : get_residue(resno)->get_phi())      // TODO: Persist previous best candidates.
+                        ((ires & 0x1) ? aa->get_psi() : aa->get_phi())      // TODO: Persist previous best candidates.
                         + angle * working_candidate[ires];
 
                     cout << candidates[ikept][icand][ires]*fiftyseven << " ";
 
-                    // TODO: Rotate and score the backbone.
-                    // Score each candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle.
+                    // Rotate the backbone.
+                    MovabilityType was_mov = aa->movability;
+                    aa->movability = MOV_ALL;
+                    if (ires & 0x1) rotate_backbone_partial(startres, endres, (dir>0) ? CA_asc : C_desc, candidates[ikept][icand][ires] - aa->get_psi());
+                    else rotate_backbone_partial(startres, endres, (dir>0) ? N_asc : CA_desc, candidates[ikept][icand][ires] - aa->get_phi());
+                    aa->movability = was_mov;
                 }
+
+                // TODO: Score the candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle.
+                anomaly = is_at[0]->get_location().get_3d_distance(should_be_at[0])
+                    + reconnect_angle_importance * is_at[1]->get_location().get_3d_distance(should_be_at[1])
+                    + reconnect_angle_importance * is_at[2]->get_location().get_3d_distance(should_be_at[2])
+                    ;
+                
+                cout << anomaly;
+
                 cout << endl;
 
                 // TODO: Keep the best reconnect_keepbest from each generation and retry.
@@ -1873,7 +1928,7 @@ float Protein::reconnect(int startres, int endres)
 
     // TODO: Use the best candidate and trial-and-error it onto connres.
 
-    // TODO: Return the anomaly.
+    return anomaly;
 }
 
 void Protein::make_helix(int startres, int endres, float phi, float psi)

From a659e216d7294d084552042169c7722e0078b8c0 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sun, 12 Nov 2023 08:11:09 -0700
Subject: [PATCH 05/17] Best candidates from each generation.

---
 src/classes/protein.cpp | 34 +++++++++++++++++++++++++++++-----
 1 file changed, 29 insertions(+), 5 deletions(-)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index 5c3b8f9f5..39e07a2df 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1849,7 +1849,7 @@ float Protein::reconnect(int startres, int endres)
     float best_score[reconnect_keepbest];
     float anomaly;
 
-    int ikept, icand, ires, jres;
+    int ikept, icand, ires, i, j, l;
     int working_candidate[length+1];
 
     Point should_be_at[4];
@@ -1873,6 +1873,12 @@ float Protein::reconnect(int startres, int endres)
 
     for (gen=0; gen<reconnect_generations; gen++)
     {
+        cout << "gen " << gen << ": ";
+        for (ikept=0; ikept < reconnect_keepbest; ikept++)
+        {
+            best_score[ikept] = -1;
+        }
+
         // TODO: try every possible combination of starting phi,psi angles +/- the angle variable.
         for (ikept=0; ikept < (gen ? reconnect_keepbest : 1); ikept++)
         {
@@ -1895,8 +1901,12 @@ float Protein::reconnect(int startres, int endres)
                         cout << "Attempt to reconnect discontinuous strand." << endl;
                         throw 0xb5712a9d;
                     }
-                    candidates[ikept][icand][ires] = 
-                        ((ires & 0x1) ? aa->get_psi() : aa->get_phi())      // TODO: Persist previous best candidates.
+                    candidates[ikept][icand][ires] =
+                        (
+                            gen
+                            ? best_candidates[ikept][ires]                      // Persist previous best candidates.
+                            : ((ires & 0x1) ? aa->get_psi() : aa->get_phi())
+                        )
                         + angle * working_candidate[ires];
 
                     cout << candidates[ikept][icand][ires]*fiftyseven << " ";
@@ -1917,9 +1927,23 @@ float Protein::reconnect(int startres, int endres)
                 
                 cout << anomaly;
 
-                cout << endl;
+                // Keep the best reconnect_keepbest from each generation.
+                for (i=0; i < reconnect_keepbest; i++)
+                {
+                    if (best_score[i] < 0 || best_score[i] > anomaly)
+                    {
+                        for (j=reconnect_keepbest-2; j>=i; j--)
+                        {
+                            best_score[j+1] = best_score[j];
+                            for (l=0; l<length; l++) best_candidates[j+1][l] = best_candidates[j][l];
+                        }
+
+                        best_score[i] = anomaly;
+                        for (l=0; l<length; l++) best_candidates[j][l] = candidates[ikept][icand][l];
+                    }
+                }
 
-                // TODO: Keep the best reconnect_keepbest from each generation and retry.
+                cout << endl;
             }
         }
 

From 352f5b77423097ea261a261f61964a404d82d0c5 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sun, 12 Nov 2023 08:32:36 -0700
Subject: [PATCH 06/17] Length bug fixes.

---
 src/classes/protein.cpp | 29 +++++++++++++++++++----------
 test/reconnect.pepd     |  2 +-
 2 files changed, 20 insertions(+), 11 deletions(-)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index 39e07a2df..f3d69bc98 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1817,8 +1817,8 @@ void Protein::conform_backbone(int startres, int endres,
 
 float Protein::reconnect(int startres, int endres)
 {
-    int length = abs(endres - startres) + 1;
-    if (length > reconnect_length_limit)
+    int strandlength = abs(endres - startres) + 1, length = strandlength*2;
+    if (strandlength > reconnect_length_limit)
     {
         cout << "Attempt to reconnect strand in excess of length limit." << endl;
         throw 0xb5712a9d;
@@ -1841,6 +1841,7 @@ float Protein::reconnect(int startres, int endres)
         throw 0xb5712a9d;
     }
 
+    char working_candidate[length+1];
     float angle = triangular;
     int gen;
     int max_candidates = pow(3, length);
@@ -1850,7 +1851,6 @@ float Protein::reconnect(int startres, int endres)
     float anomaly;
 
     int ikept, icand, ires, i, j, l;
-    int working_candidate[length+1];
 
     Point should_be_at[4];
     Atom* is_at[4];
@@ -1873,7 +1873,6 @@ float Protein::reconnect(int startres, int endres)
 
     for (gen=0; gen<reconnect_generations; gen++)
     {
-        cout << "gen " << gen << ": ";
         for (ikept=0; ikept < reconnect_keepbest; ikept++)
         {
             best_score[ikept] = -1;
@@ -1885,12 +1884,7 @@ float Protein::reconnect(int startres, int endres)
             for (ires = 0; ires < length; ires++) working_candidate[ires] = -1;
             for (icand = 0; icand < max_candidates; icand++)
             {
-                for (ires = 0; ires < length; ires++)
-                {
-                    working_candidate[ires]++;
-                    if (working_candidate[ires] > 1) working_candidate[ires] = -1;
-                    else break;
-                }
+                cout << "gen " << gen << " kept " << ikept << " iter " << icand << " ";
 
                 for (ires = 0; ires < length; ires++)
                 {
@@ -1901,6 +1895,7 @@ float Protein::reconnect(int startres, int endres)
                         cout << "Attempt to reconnect discontinuous strand." << endl;
                         throw 0xb5712a9d;
                     }
+
                     candidates[ikept][icand][ires] =
                         (
                             gen
@@ -1942,11 +1937,25 @@ float Protein::reconnect(int startres, int endres)
                         for (l=0; l<length; l++) best_candidates[j][l] = candidates[ikept][icand][l];
                     }
                 }
+                
+                for (ires = 0; ires < length; ires++)
+                {
+                    working_candidate[ires]++;
+                    if (working_candidate[ires] > 1) working_candidate[ires] = -1;
+                    else break;
+                }
 
                 cout << endl;
             }
         }
 
+        cout << endl << "Best from generation:" << endl;
+        for (i=0; i < reconnect_keepbest; i++)
+        {
+            cout << i+1 << ": " << best_score[i] << endl;
+        }
+        cout << endl;
+
         angle /= 2;
     }
 
diff --git a/test/reconnect.pepd b/test/reconnect.pepd
index d99fe6a7c..453750662 100644
--- a/test/reconnect.pepd
+++ b/test/reconnect.pepd
@@ -5,4 +5,4 @@ LET @axis = @5.50 - @6.48
 ECHO @axis
 
 ROTATE @axis -15 %6.48 %6.28 %6.59
-CONNECT %6.63 %6.59
+CONNECT %6.61 %6.59

From ae8133fdcc81f0a4d5ae8f318605f09659310042 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sun, 12 Nov 2023 09:05:35 -0700
Subject: [PATCH 07/17] Bug fix for keep best candidates.

---
 src/classes/protein.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index f3d69bc98..eee0024ff 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1935,6 +1935,8 @@ float Protein::reconnect(int startres, int endres)
 
                         best_score[i] = anomaly;
                         for (l=0; l<length; l++) best_candidates[j][l] = candidates[ikept][icand][l];
+
+                        break;
                     }
                 }
                 

From ee0f67607f283cfc709fb47ee1c0e83bc9814506 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sun, 12 Nov 2023 09:37:35 -0700
Subject: [PATCH 08/17] Debug show angles.

---
 src/classes/protein.cpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index eee0024ff..bb938ef0e 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1954,7 +1954,9 @@ float Protein::reconnect(int startres, int endres)
         cout << endl << "Best from generation:" << endl;
         for (i=0; i < reconnect_keepbest; i++)
         {
-            cout << i+1 << ": " << best_score[i] << endl;
+            cout << i+1 << ": ";
+            for (l=0; l<length; l++) cout << best_candidates[i][l]*fiftyseven << " ";
+            cout << "= " << best_score[i] << endl;
         }
         cout << endl;
 

From 9aafa61761eaf32bd9c16f999b2f93ed0890d25e Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sun, 12 Nov 2023 12:24:36 -0700
Subject: [PATCH 09/17] Better chance of success.

---
 test/reconnect.pepd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/reconnect.pepd b/test/reconnect.pepd
index 453750662..82d14010a 100644
--- a/test/reconnect.pepd
+++ b/test/reconnect.pepd
@@ -4,5 +4,5 @@ LOAD "pdbs/OR51/OR51E2.upright.pdb"
 LET @axis = @5.50 - @6.48
 ECHO @axis
 
-ROTATE @axis -15 %6.48 %6.28 %6.59
+ROTATE @axis -5 %6.48 %6.28 %6.59
 CONNECT %6.61 %6.59

From 405a63a33c53dac80c4ed0e467be926e2ab9b429 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Mon, 13 Nov 2023 22:35:04 -0700
Subject: [PATCH 10/17] Now it makes a much faster mistake.

---
 src/classes/protein.cpp | 77 +++++++++++++++++++++++++++++++++++++----
 src/classes/protein.h   |  2 +-
 test/reconnect.pepd     |  3 ++
 3 files changed, 74 insertions(+), 8 deletions(-)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index bb938ef0e..d2cf1f197 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1554,7 +1554,7 @@ void Protein::rotate_backbone(int resno, bb_rot_dir dir, float angle)
     }
 }
 
-void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir, float angle)
+void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir, float angle, bool bbao)
 {
     save_undo_state();
     if (startres == endres) return;
@@ -1568,8 +1568,50 @@ void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir,
 
     AminoAcid* bendy = get_residue(startres);
     if (!bendy) return;
-    LocatedVector lv = bendy->rotate_backbone(dir, angle);
-    bendy->ensure_pi_atoms_coplanar();
+    LocatedVector lv;
+    if (bbao)
+    {
+        Atom *N = bendy->get_atom("N"), *H = bendy->HN_or_substitute(), *CA = bendy->get_atom("CA"), *C = bendy->get_atom("C"), *O = bendy->get_atom("O");
+        Bond* b;
+        switch(dir)
+        {
+            case N_asc:
+            b = N->get_bond_between(CA);
+            lv = (SCoord)CA->get_location().subtract(N->get_location());
+            lv.origin = N->get_location();
+            b->rotate(angle, true, true);
+            break;
+
+            case CA_desc:
+            b = CA->get_bond_between(N);
+            lv = (SCoord)CA->get_location().subtract(N->get_location());
+            lv.origin = CA->get_location();
+            b->rotate(angle, true, true);
+            break;
+
+            case CA_asc:
+            b = CA->get_bond_between(C);
+            lv = (SCoord)CA->get_location().subtract(C->get_location());
+            lv.origin = CA->get_location();
+            b->rotate(angle, true, true);
+            break;
+
+            case C_desc:
+            b = C->get_bond_between(CA);
+            lv = (SCoord)C->get_location().subtract(CA->get_location());
+            lv.origin = C->get_location();
+            b->rotate(angle, true, true);
+            break;
+
+            default:
+            ;
+        }
+    }
+    else
+    {
+        lv = bendy->rotate_backbone(dir, angle);
+        bendy->ensure_pi_atoms_coplanar();
+    }
 
     if (lv.r)
     {
@@ -1579,12 +1621,12 @@ void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir,
         for (i=startres+inc; movable = get_residue(i); i+=inc)
         {
             movable->rotate(lv, angle);
-            movable->ensure_pi_atoms_coplanar();
+            if (!bbao) movable->ensure_pi_atoms_coplanar();
             if (i == endres) break;
         }
     }
 
-    set_clashables();
+    if (!bbao) set_clashables();
 }
 
 void Protein::conform_backbone(int startres, int endres, int iters, bool backbone_atoms_only)
@@ -1909,8 +1951,12 @@ float Protein::reconnect(int startres, int endres)
                     // Rotate the backbone.
                     MovabilityType was_mov = aa->movability;
                     aa->movability = MOV_ALL;
-                    if (ires & 0x1) rotate_backbone_partial(startres, endres, (dir>0) ? CA_asc : C_desc, candidates[ikept][icand][ires] - aa->get_psi());
-                    else rotate_backbone_partial(startres, endres, (dir>0) ? N_asc : CA_desc, candidates[ikept][icand][ires] - aa->get_phi());
+                    if (ires & 0x1) rotate_backbone_partial(resno, endres,
+                        (dir>0) ? CA_asc : C_desc,
+                        candidates[ikept][icand][ires] - aa->get_psi(), true);
+                    else rotate_backbone_partial(resno, endres,
+                        (dir>0) ? N_asc : CA_desc,
+                        candidates[ikept][icand][ires] - aa->get_phi(), true);
                     aa->movability = was_mov;
                 }
 
@@ -1964,7 +2010,24 @@ float Protein::reconnect(int startres, int endres)
     }
 
     // TODO: Use the best candidate and trial-and-error it onto connres.
+    for (ires=0; ires<length; ires++)
+    {
+        int resno = startres + dir*(ires>>1);
+        AminoAcid* aa = get_residue(resno);
 
+        // Rotate the backbone.
+        MovabilityType was_mov = aa->movability;
+        aa->movability = MOV_ALL;
+        if (ires & 0x1) rotate_backbone_partial(resno, endres,
+            (dir>0) ? CA_asc : C_desc,
+            best_candidates[0][ires] - aa->get_psi(), true);
+        else rotate_backbone_partial(resno, endres,
+            (dir>0) ? N_asc : CA_desc,
+            best_candidates[0][ires] - aa->get_phi(), true);
+        aa->movability = was_mov;
+    }
+
+    set_clashables();
     return anomaly;
 }
 
diff --git a/src/classes/protein.h b/src/classes/protein.h
index 5f529b335..58f26d990 100755
--- a/src/classes/protein.h
+++ b/src/classes/protein.h
@@ -143,8 +143,8 @@ class Protein
     LocRotation rotate_piece(int start_res, int end_res, Point origin, SCoord axis, float theta);
 
     void rotate_backbone(int residue_no, bb_rot_dir direction, float angle);
+    void rotate_backbone_partial(int startres, int endres, bb_rot_dir direction, float angle, bool backbone_atoms_only = false);
     void conform_backbone(int startres, int endres, Atom* a, Point target, int iters = 50);
-    void rotate_backbone_partial(int startres, int endres, bb_rot_dir direction, float angle);
     void conform_backbone(int startres, int endres, int iters = 50, bool backbone_atoms_only = false);
     void conform_backbone(int startres, int endres, Atom* a1, Point target1, Atom* a2, Point target2, Atom* a3, Point target3, int iters = 50);
     void conform_backbone(int startres, int endres, Atom* a1, Point target1, Atom* a2, Point target2, int iters = 50, bool backbone_atoms_only = false);
diff --git a/test/reconnect.pepd b/test/reconnect.pepd
index 82d14010a..e4f324bbd 100644
--- a/test/reconnect.pepd
+++ b/test/reconnect.pepd
@@ -6,3 +6,6 @@ ECHO @axis
 
 ROTATE @axis -5 %6.48 %6.28 %6.59
 CONNECT %6.61 %6.59
+
+SAVE "tmp/OR51E2.reconnected.pdb"
+

From 1d19eb8dcbd3d35a763924dca7969619584b62eb Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Wed, 15 Nov 2023 10:07:42 -0700
Subject: [PATCH 11/17] Constants.

---
 src/classes/constants.h | 4 ++--
 src/classes/point.cpp   | 2 +-
 test/reconnect.pepd     | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/classes/constants.h b/src/classes/constants.h
index 622a8d2a6..b7908c7e7 100755
--- a/src/classes/constants.h
+++ b/src/classes/constants.h
@@ -319,8 +319,8 @@
 #define homology_region_optimization 0
 
 #define reconnect_length_limit 5
-#define reconnect_generations 8
-#define reconnect_keepbest 4
+#define reconnect_generations 10
+#define reconnect_keepbest 5
 #define reconnect_angle_importance 0.1
 
 //////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/src/classes/point.cpp b/src/classes/point.cpp
index ad766d631..b341e0d69 100755
--- a/src/classes/point.cpp
+++ b/src/classes/point.cpp
@@ -301,7 +301,7 @@ float find_angle_along_vector(Point* pt1, Point* pt2, Point* source, SCoord* v)
     Point lpt2 = pt2->subtract(source);
 
     // Rotate points so v becomes Z axis.
-    Point cen;
+    Point cen(0,0,0);
     Point z(0,0,1);
     Rotation rots = align_points_3d(&vp, &z, &cen);
     Point npt1 = rotate3D(&lpt1, &cen, &rots.v, rots.a);
diff --git a/test/reconnect.pepd b/test/reconnect.pepd
index e4f324bbd..f6068f7b0 100644
--- a/test/reconnect.pepd
+++ b/test/reconnect.pepd
@@ -5,7 +5,7 @@ LET @axis = @5.50 - @6.48
 ECHO @axis
 
 ROTATE @axis -5 %6.48 %6.28 %6.59
-CONNECT %6.61 %6.59
+CONNECT %6.63 %6.59
 
 SAVE "tmp/OR51E2.reconnected.pdb"
 

From ce5f4839666b6d79d58c485e28a96b9c4f02a128 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Wed, 15 Nov 2023 20:17:04 -0700
Subject: [PATCH 12/17] Helix option in Ramachandran plot.

---
 src/ramachandran.cpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/ramachandran.cpp b/src/ramachandran.cpp
index afd04932b..d7a9b2578 100644
--- a/src/ramachandran.cpp
+++ b/src/ramachandran.cpp
@@ -9,11 +9,13 @@ int main(int argc, char** argv)
     int i;
     FILE* fp;
     bool show_colors = true;
+    bool skip_helices = false;
     for (i=1; i<argc; i++)
     {
         if (argv[i][0] == '-')
         {
             if (!strcmp(argv[i], "-n") || !strcmp(argv[i], "--numbers")) show_colors = false;
+            if (!strcmp(argv[i], "-h") || !strcmp(argv[i], "--skip-helices")) skip_helices = true;
             continue;
         }
 
@@ -37,6 +39,7 @@ int main(int argc, char** argv)
     {
         AminoAcid* aa = p.get_residue(i);
         if (!aa) continue;
+        if (skip_helices && aa->is_alpha_helix()) continue;
         float phi = aa->get_phi(), psi = aa->get_psi();
 
         phi = fmod(phi, (M_PI*2)); if (phi >= M_PI) phi -= M_PI*2;

From fd237245a45516dcf0073328885f208f13cf0593 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Wed, 15 Nov 2023 22:13:04 -0700
Subject: [PATCH 13/17] More Ramachandran options.

---
 src/ramachandran.cpp | 34 ++++++++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/src/ramachandran.cpp b/src/ramachandran.cpp
index d7a9b2578..916b7657a 100644
--- a/src/ramachandran.cpp
+++ b/src/ramachandran.cpp
@@ -10,12 +10,22 @@ int main(int argc, char** argv)
     FILE* fp;
     bool show_colors = true;
     bool skip_helices = false;
+    bool include_glycine = true, include_not_glycine = true;
+    bool include_proline = true, include_not_proline = true;
+    bool include_pre_proline = true, include_not_pre_proline = true;
     for (i=1; i<argc; i++)
     {
         if (argv[i][0] == '-')
         {
             if (!strcmp(argv[i], "-n") || !strcmp(argv[i], "--numbers")) show_colors = false;
             if (!strcmp(argv[i], "-h") || !strcmp(argv[i], "--skip-helices")) skip_helices = true;
+            if (!strcmp(argv[i], "-g") || !strcmp(argv[i], "--omit-glycine")) include_glycine = false;
+            if (!strcmp(argv[i], "-G") || !strcmp(argv[i], "--only-glycine")) include_not_glycine = false;
+            if (!strcmp(argv[i], "-p") || !strcmp(argv[i], "--omit-proline")) include_proline = false;
+            if (!strcmp(argv[i], "-P") || !strcmp(argv[i], "--only-proline")) include_not_proline = false;
+            if (!strcmp(argv[i], "-b") || !strcmp(argv[i], "--omit-pre-proline")) include_pre_proline = false;
+            if (!strcmp(argv[i], "-B") || !strcmp(argv[i], "--only-pre-proline")) include_not_pre_proline = false;
+            if (!include_glycine && !include_not_glycine) return -1;
             continue;
         }
 
@@ -40,6 +50,30 @@ int main(int argc, char** argv)
         AminoAcid* aa = p.get_residue(i);
         if (!aa) continue;
         if (skip_helices && aa->is_alpha_helix()) continue;
+        if (aa->get_letter() == 'G')
+        {
+            if (!include_glycine || !include_not_proline) continue;
+        }
+        else if (aa->get_letter() == 'P')
+        {
+            if (!include_proline || !include_not_glycine) continue;
+        }
+        else if (!include_not_glycine || !include_not_proline) continue;
+
+        if (!include_pre_proline || !include_not_pre_proline)
+        {
+            AminoAcid* next = p.get_residue(i+1);
+            if (next)
+            {
+                if (next->get_letter() == 'P')
+                {
+                    if (!include_pre_proline) continue;
+                }
+                else if (!include_not_pre_proline) continue;
+            }
+        }
+
+
         float phi = aa->get_phi(), psi = aa->get_psi();
 
         phi = fmod(phi, (M_PI*2)); if (phi >= M_PI) phi -= M_PI*2;

From a84f8b37a67e4a938873513c75f1c0ddbeddcea1 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Wed, 15 Nov 2023 23:08:16 -0700
Subject: [PATCH 14/17] Directionality.

---
 src/classes/protein.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index d2cf1f197..a123198ca 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -1953,10 +1953,10 @@ float Protein::reconnect(int startres, int endres)
                     aa->movability = MOV_ALL;
                     if (ires & 0x1) rotate_backbone_partial(resno, endres,
                         (dir>0) ? CA_asc : C_desc,
-                        candidates[ikept][icand][ires] - aa->get_psi(), true);
+                        (candidates[ikept][icand][ires] - aa->get_psi()) * dir, true);
                     else rotate_backbone_partial(resno, endres,
                         (dir>0) ? N_asc : CA_desc,
-                        candidates[ikept][icand][ires] - aa->get_phi(), true);
+                        (candidates[ikept][icand][ires] - aa->get_phi()) * dir, true);
                     aa->movability = was_mov;
                 }
 
@@ -2020,10 +2020,10 @@ float Protein::reconnect(int startres, int endres)
         aa->movability = MOV_ALL;
         if (ires & 0x1) rotate_backbone_partial(resno, endres,
             (dir>0) ? CA_asc : C_desc,
-            best_candidates[0][ires] - aa->get_psi(), true);
+            (best_candidates[0][ires] - aa->get_psi()) * dir, true);
         else rotate_backbone_partial(resno, endres,
             (dir>0) ? N_asc : CA_desc,
-            best_candidates[0][ires] - aa->get_phi(), true);
+            (best_candidates[0][ires] - aa->get_phi()) * dir, true);
         aa->movability = was_mov;
     }
 

From 7d16cd6164a205b7e777b722a998b63a539e1756 Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Thu, 16 Nov 2023 17:29:20 -0700
Subject: [PATCH 15/17] ScorePDB fixes.

---
 src/classes/atom.cpp     |  2 +-
 src/classes/molecule.cpp | 33 ++++++++++++++++++++++++++++++---
 src/classes/scoring.cpp  |  2 +-
 src/score_pdb.cpp        |  7 ++++---
 4 files changed, 36 insertions(+), 8 deletions(-)

diff --git a/src/classes/atom.cpp b/src/classes/atom.cpp
index 59f6ce583..448304927 100755
--- a/src/classes/atom.cpp
+++ b/src/classes/atom.cpp
@@ -572,7 +572,7 @@ float Atom::get_charge()
 {
     if (Z == 1)
     {
-        if (bonded_to && bonded_to[0].btom)
+        if (bonded_to[0].btom)
         {
             float bchg = bonded_to[0].btom->charge;
             if (!bchg) bchg = bonded_to[0].btom->is_conjugated_to_charge();
diff --git a/src/classes/molecule.cpp b/src/classes/molecule.cpp
index a9622f7c3..65927a156 100644
--- a/src/classes/molecule.cpp
+++ b/src/classes/molecule.cpp
@@ -423,6 +423,7 @@ void Molecule::hydrogenate(bool steric_only)
 
         float valence = atoms[i]->get_valence();
         if (valence > 4) valence = 8 - valence;
+        if (valence > 1) valence += atoms[i]->get_charge();
 
         // cout << atoms[i]->name << " has valence " << valence << endl;
 
@@ -461,7 +462,7 @@ void Molecule::hydrogenate(bool steric_only)
             char hname[5];
             sprintf(hname, "H%d", atcount+1);
             Atom* H = add_atom("H", hname, atoms[i], 1);
-            //cout << "Adding " << hname << " to " << atoms[i]->name << " whose valence is " << valence << " and has " << bcardsum << " bonds already." << endl;
+            // if (!is_residue()) cout << "Adding " << hname << " to " << atoms[i]->name << " whose valence is " << valence << " and has " << bcardsum << " bonds already." << endl;
 
             /*atoms[i]->clear_geometry_cache();
             SCoord v = atoms[i]->get_next_free_geometry(1);
@@ -879,10 +880,16 @@ int Molecule::from_pdb(FILE* is, bool het_only)
                     }
                     esym[1] &= 0x5f;
 
-                    Point aloc(atof(words[5]), atof(words[6]),atof(words[7]));
+                    int offset = 0;
+                    if (words[4][0] >= 'A' || words[5] < &buffer[29]) offset++;
+
+                    Point aloc(atof(words[5+offset]), atof(words[6+offset]),atof(words[7+offset]));
+
+                    // cout << esym << " " << aloc << endl;
 
                     Atom* a = add_atom(esym, words[2], &aloc, 0, 0, charge);
                     added++;
+                    a->pdbidx = atoi(words[1]);
 
                     // a->residue = atoi(words[4]);
 
@@ -903,6 +910,20 @@ int Molecule::from_pdb(FILE* is, bool het_only)
                     ;
                 }
             }
+            else if (!strcmp(words[0], "CONECT"))
+            {
+                int i, i1 = atoi(words[1]), i2 = atoi(words[2]);
+                Atom *a = nullptr, *b = nullptr;
+
+                for (i=0; atoms[i]; i++)
+                {
+                    if (atoms[i]->pdbidx == i1) a = atoms[i];
+                    if (atoms[i]->pdbidx == i2) b = atoms[i];
+                    if (a && b) break;
+                }
+
+                if (a && b) a->bond_to(b, 1);              // TODO: Detect coplanar atoms and add pi bonds where found.
+            }
         }
         buffer[0] = 0;
 
@@ -2402,7 +2423,7 @@ float Molecule::get_intermol_binding(Molecule** ligands, bool subtract_clashes)
                     ligands[l]->atcount = j;
                     break;
                 }
-                if (atoms[i]->is_backbone && ligands[l]->atoms[j]->is_backbone
+                if (atoms[i]->is_backbone && atoms[i]->residue && ligands[l]->atoms[j]->is_backbone && ligands[l]->atoms[j]->residue
                     &&
                     (	(	atoms[i]->residue == ligands[l]->atoms[j]->residue - 1
                             &&
@@ -2421,14 +2442,20 @@ float Molecule::get_intermol_binding(Molecule** ligands, bool subtract_clashes)
                 float r = ligands[l]->atoms[j]->get_location().get_3d_distance(&aloc);
                 if (r < _INTERA_R_CUTOFF)
                 {
+                    // cout << "🐝 ";
                     if (	!shielded(atoms[i], ligands[l]->atoms[j])
                             &&
                             !ligands[l]->shielded(atoms[i], ligands[l]->atoms[j])
                        )
                     {
                         missed_connection.r = 0;
+                        // if (ligands[l]->atoms[j]->residue) cout << ligands[l]->atoms[j]->residue << ":";
+                        // cout << ligands[l]->atoms[j]->name << " " << ligands[l]->atoms[j]->get_charge() << " ~ ";
+                        // if (atoms[i]->residue) cout << atoms[i]->residue << ":";
+                        // cout << atoms[i]->name << " " << atoms[i]->get_charge() << " ";
                         // cout << ligands[l]->atoms[j]->get_location().subtract(atoms[i]->get_location()) << ": ";
                         float abind = InteratomicForce::total_binding(atoms[i], ligands[l]->atoms[j]);
+                        // cout << abind << endl;
                         if (abind && !isnan(abind) && !isinf(abind))
                         {
                             if (abind > 0 && minimum_searching_aniso && ligands[l]->priority) abind *= 1.5;
diff --git a/src/classes/scoring.cpp b/src/classes/scoring.cpp
index 909b8bd02..549a5e590 100644
--- a/src/classes/scoring.cpp
+++ b/src/classes/scoring.cpp
@@ -228,7 +228,7 @@ DockResult::DockResult(Protein* protein, Molecule* ligand, Point size, int* addl
         #endif
 
         sprintf(metrics[metcount], "%s%d", reaches_spheroid[i]->get_3letter(), resno);
-        // cout << metrics[metcount] << ": " << lb << " . ";
+        // cout << metrics[metcount] << ": " << lb << ". ";
 
         if (differential_dock)
         {
diff --git a/src/score_pdb.cpp b/src/score_pdb.cpp
index 8451dd31d..a52929756 100644
--- a/src/score_pdb.cpp
+++ b/src/score_pdb.cpp
@@ -10,6 +10,8 @@ using namespace std;
 
 int main(int argc, char** argv)
 {
+    _INTERA_R_CUTOFF = _DEFAULT_INTERA_R_CUTOFF;
+
     if (argc < 2)
     {
         cout << "Usage:" << endl << endl;
@@ -34,13 +36,12 @@ int main(int argc, char** argv)
     fseek(fp, 0, 0);
     m.from_pdb(fp, true);
 
-    AminoAcid* aa = p.get_residue_bw(3, 50);
-    if (!aa) aa = p.get_residue_bw(6, 48);
-    if (!aa) aa = p.get_residue(100);
+    AminoAcid* aa = p.get_residue(100);
     if (!aa) aa = p.get_residue(p.get_end_resno()/2);
 
     if (aa && aa->get_hydrogen_count() < 3)
     {
+        cout << "Hydrogenating..." << endl;
         int i, n = p.get_end_resno();
         for (i=0; i<n; i++)
         {

From 91253dad38e8bb2f8f5b3b266dda3db4c0aaf2ca Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Fri, 17 Nov 2023 20:12:47 -0700
Subject: [PATCH 16/17] BWCOPY command.

---
 SCRIPTING.md               | 13 ++++++++++++-
 data/taar_activations.pepd | 17 +++++++++++++++++
 src/classes/constants.h    |  1 +
 src/classes/molecule.cpp   | 16 +++++++---------
 src/classes/protein.cpp    |  7 +++++++
 src/classes/protein.h      |  1 +
 src/interpreter.cpp        | 23 ++++++++++++++++++++++-
 7 files changed, 67 insertions(+), 11 deletions(-)
 create mode 100644 data/taar_activations.pepd

diff --git a/SCRIPTING.md b/SCRIPTING.md
index 2d55526b3..6dac315b2 100644
--- a/SCRIPTING.md
+++ b/SCRIPTING.md
@@ -181,11 +181,22 @@ BWCENTER
 ```
 
 If the current working strand is a seven-helix protein (7HP), with Ballesteros-Weinstein numbering for all its transmembrane helices,
-then `BWCENTER` will obtain the center of all n.50:CA atom locations and recenter the entire protein.
+then `BWCENTER` will obtain the center of all n.50:CA atom locations and recenter the entire strand.
 This is useful for comparing various active and inactive states of GPCRs; the 1.50 through 7.50 residues will be assumed to be the most stationary
 parts of the protein, and variations in structure can be observed with minimal global transformational anomalies.
 
 
+# BWCOPY
+Example:
+```
+BWCOPY A B
+```
+
+Copies the transmembrane helix regions and Ballesteros-Weinstein numbers from one strand to another. This is useful if working with two models of
+the same protein, one of which is missing the region definitions, for example if one of the models is included with PrimaryDock and the other is not.
+The first parameter is the source strand; the second is the destination.
+
+
 # CANMOVE
 Example:
 ```
diff --git a/data/taar_activations.pepd b/data/taar_activations.pepd
new file mode 100644
index 000000000..cc400f09e
--- /dev/null
+++ b/data/taar_activations.pepd
@@ -0,0 +1,17 @@
+
+LOAD "pdbs/TAAR/TAAR1.upright.pdb" A I
+CENTER
+UPRIGHT I
+BWCENTER
+
+LOAD "pdbs/8jlj.pdb" R A
+CENTER
+
+BWCOPY I A
+UPRIGHT A
+BWCENTER
+
+# DUMP QUIT
+
+SAVE "pdbs/TAAR1.3-iodothyronamine.pdb"
+
diff --git a/src/classes/constants.h b/src/classes/constants.h
index b7908c7e7..90fd9fce5 100755
--- a/src/classes/constants.h
+++ b/src/classes/constants.h
@@ -70,6 +70,7 @@
 #define _ALLOW_FLEX_RINGS 0
 #define _shield_angle (130.0 * fiftyseventh)
 #define _shield_angle_pi (100.0 * fiftyseventh)
+#define _shield_angle_opposite_charge (65.0 * fiftyseventh)
 #define _can_clash_angle (180.0 * fiftyseventh)
 #define _fullrot_stepdeg 20
 #define _fullrot_steprad (fiftyseventh*_fullrot_stepdeg)
diff --git a/src/classes/molecule.cpp b/src/classes/molecule.cpp
index 65927a156..7e4c2d052 100644
--- a/src/classes/molecule.cpp
+++ b/src/classes/molecule.cpp
@@ -2213,9 +2213,11 @@ bool Molecule::shielded(Atom* a, Atom* b) const
     a->shielding_angle = b->shielding_angle = 0;
 
     Point aloc = a->get_location(), bloc = b->get_location();
+    float achg = a->get_charge(), bchg = b->get_charge();
     for (i=0; atoms[i]; i++)
     {
         Atom* ai = atoms[i];
+        float aichg = ai->get_charge();
         if (!ai) break;
         if (ai == a || ai == b) continue;
         float rai = ai->distance_to(a);
@@ -2228,15 +2230,11 @@ bool Molecule::shielded(Atom* a, Atom* b) const
         if (f3da > a->shielding_angle) a->shielding_angle = b->shielding_angle = f3da;
         if (f3da > _shield_angle)
         {
-            if (last_iter && (a->residue == 114 || b->residue == 114) && ((a->residue + b->residue) == 114))
-            {
-                /*cout << ai->name << " shields "
-                	 << a->residue << ":" << a->name << "..."
-                	 << b->residue << ":" << b->name
-                	 << " angle " << (f3da*fiftyseven)
-                	 << endl;*/
-                return true;
-            }
+            return true;
+        }
+        if (sgn(achg) == sgn(bchg) && sgn(achg) == -sgn(aichg) && f3da > _shield_angle_opposite_charge)
+        {
+            return true;
         }
     }
 
diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp
index a123198ca..56e6d9799 100755
--- a/src/classes/protein.cpp
+++ b/src/classes/protein.cpp
@@ -993,6 +993,13 @@ int Protein::get_bw50(int helixno)
     return Ballesteros_Weinstein[helixno];
 }
 
+void Protein::set_bw50(int hxno, int resno)
+{
+    if (hxno < 1 || hxno > 78) return;
+    if (Ballesteros_Weinstein[hxno]) return;
+    Ballesteros_Weinstein[hxno] = resno;
+}
+
 int Protein::get_seq_length()
 {
     if (!sequence) return 0;
diff --git a/src/classes/protein.h b/src/classes/protein.h
index 58f26d990..5d90b155c 100755
--- a/src/classes/protein.h
+++ b/src/classes/protein.h
@@ -61,6 +61,7 @@ class Protein
     void delete_sidechains(int startres, int endres);
     MetalCoord* coordinate_metal(Atom* metal, int residues, int* resnos, std::vector<string> res_anames);
     void set_region(std::string name, int start, int end);
+    void set_bw50(int helixno, int resno);
     void renumber_residues(int startres, int endres, int new_startres);
     bool disulfide_bond(int resno1, int resno2);
 
diff --git a/src/interpreter.cpp b/src/interpreter.cpp
index 6cdc012fd..7e8192990 100644
--- a/src/interpreter.cpp
+++ b/src/interpreter.cpp
@@ -964,7 +964,28 @@ int main(int argc, char** argv)
                 SCoord new_center = old_center;
                 new_center.r *= -1;
                 working->move_piece(1, working->get_end_resno(), new_center);
-            }   // BWCENTER
+            }	// BWCENTER
+
+            else if (!strcmp(words[0], "BWCOPY"))
+            {
+                if (!words[1] || !words[2]) raise_error("Insufficient parameters given for BWCOPY.");
+
+                char c = words[1][0];
+                if (c < 'A' || c > 'Z') raise_error("Invalid source strand given for BWCOPY.");
+                Protein* source = strands[c-'A'];
+
+                c = words[2][0];
+                if (c < 'A' || c > 'Z') raise_error("Invalid destination strand given for BWCOPY.");
+                Protein* dest = strands[c-'A'];
+
+                for (l=1; l<=8; l++)
+                {
+                    std::string rgname = (std::string)"TMR" + std::to_string(l);
+                    int sr = source->get_region_start(rgname), er = source->get_region_end(rgname), bw50 = source->get_bw50(l);
+                    if (sr > 0 && er > 0) dest->set_region(rgname, sr, er);
+                    if (bw50 > 0) dest->set_bw50(l, bw50);
+                }
+            }   // BWCOPY
 
             else if (!strcmp(words[0], "CANMOVE"))
             {

From 4be27c83e5345c1fdb460ac3749c93fdf29f834c Mon Sep 17 00:00:00 2001
From: "PrimaryOdors.org" <julie@primaryodors.org>
Date: Sat, 18 Nov 2023 10:25:53 -0700
Subject: [PATCH 17/17] Direction of test.

---
 data/taar_activations.pepd | 44 +++++++++++++++++++++++++++++++++-----
 test/reconnect.pepd        |  3 ++-
 2 files changed, 41 insertions(+), 6 deletions(-)

diff --git a/data/taar_activations.pepd b/data/taar_activations.pepd
index cc400f09e..67965306b 100644
--- a/data/taar_activations.pepd
+++ b/data/taar_activations.pepd
@@ -1,17 +1,51 @@
 
 LOAD "pdbs/TAAR/TAAR1.upright.pdb" A I
-CENTER
 UPRIGHT I
 BWCENTER
 
-LOAD "pdbs/8jlj.pdb" R A
-CENTER
+DOWNLOAD RCSB 8JLN "pdbs/8jln.pdb" ONCE
+LOAD "pdbs/8jln.pdb" R A
+BWCOPY I A
+UPRIGHT A
+BWCENTER
+SAVE "pdbs/TAAR1.3-iodothyronamine.pdb"
+
 
+DOWNLOAD RCSB 8JLO "pdbs/8jlo.pdb" ONCE
+LOAD "pdbs/8jlo.pdb" R A
 BWCOPY I A
 UPRIGHT A
 BWCENTER
+SAVE "pdbs/TAAR1.ulotaront.pdb"
 
-# DUMP QUIT
 
-SAVE "pdbs/TAAR1.3-iodothyronamine.pdb"
+DOWNLOAD RCSB 8JLP "pdbs/8jlp.pdb" ONCE
+LOAD "pdbs/8jlp.pdb" R A
+BWCOPY I A
+UPRIGHT A
+BWCENTER
+SAVE "pdbs/TAAR1.ralmitaront.pdb"
 
+
+DOWNLOAD RCSB 8JLQ "pdbs/8jlq.pdb" ONCE
+LOAD "pdbs/8jlq.pdb" R A
+BWCOPY I A
+UPRIGHT A
+BWCENTER
+SAVE "pdbs/TAAR1.fenoldopam.pdb"
+
+
+DOWNLOAD RCSB 8JLR "pdbs/8jlr.pdb" ONCE
+LOAD "pdbs/8jlr.pdb" R A
+BWCOPY I A
+UPRIGHT A
+BWCENTER
+SAVE "pdbs/TAAR1.A77636.pdb"
+
+
+DOWNLOAD RCSB 8JSO "pdbs/8jso.pdb" ONCE
+LOAD "pdbs/8jso.pdb" R A
+BWCOPY I A
+UPRIGHT A
+BWCENTER
+SAVE "pdbs/TAAR1.amphetamine.pdb"
diff --git a/test/reconnect.pepd b/test/reconnect.pepd
index f6068f7b0..a7db1db69 100644
--- a/test/reconnect.pepd
+++ b/test/reconnect.pepd
@@ -5,7 +5,8 @@ LET @axis = @5.50 - @6.48
 ECHO @axis
 
 ROTATE @axis -5 %6.48 %6.28 %6.59
-CONNECT %6.63 %6.59
+CONNECT %5.68 %6.28
+# CONNECT %6.63 %6.59
 
 SAVE "tmp/OR51E2.reconnected.pdb"