Broken strand reconnection

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:
 ```

data/taar_activations.pepd

@@ -0,0 +1,51 @@
+
+LOAD "pdbs/TAAR/TAAR1.upright.pdb" A I
+UPRIGHT I
+BWCENTER
+
+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"
+
+
+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"

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;
 }

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);

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();

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)
@@ -318,6 +319,11 @@
 #define homology_long_axis_rotations 0
 #define homology_region_optimization 0
 
+#define reconnect_length_limit 5
+#define reconnect_generations 10
+#define reconnect_keepbest 5
+#define reconnect_angle_importance 0.1
+
 //////////////////////////////////////////////////////////////////////////////////////////////
 
 // Debugging stuff.

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;
 
@@ -2192,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);
@@ -2207,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;
         }
     }
 
@@ -2402,7 +2421,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 +2440,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;

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);