From 442d7ca981553bcce5f7b7b0546f62266387dbc0 Mon Sep 17 00:00:00 2001 From: Anton Butsev Date: Mon, 19 Sep 2022 15:57:50 -0400 Subject: [PATCH] View/color salt-bridge interactions differently from hydrogen bonds (PYMOL-3846) --- data/setting_help.csv | 1 + layer1/SettingInfo.h | 1 + layer2/ObjectDist.cpp | 3 + layer3/Executive.cpp | 5 +- layer3/Interactions.cpp | 313 +++++++++++++++++++++++++++++++++------- layer3/Interactions.h | 33 +++++ modules/pymol/menu.py | 9 ++ 7 files changed, 311 insertions(+), 54 deletions(-) diff --git a/data/setting_help.csv b/data/setting_help.csv index 186b9e374..19a678fb6 100644 --- a/data/setting_help.csv +++ b/data/setting_help.csv @@ -298,6 +298,7 @@ sdf (default)","string","mol2","0" "halogen_bond_as_acceptor_min_donor_angle","Halogen-bonds as acceptor minumum donor angle.","float","120.0","0" "halogen_bond_as_acceptor_min_acceptor_angle","Halogen-bonds as acceptor minumum acceptor angle.","float","90.0","0" "halogen_bond_as_acceptor_max_acceptor_angle","Halogen-bonds as acceptor maximum acceptor angle.","float","170.0","0" +"salt_bridge_distance","Salt-bridge maximum distance.","float","5.0","0" "half_bonds","controls whether or not half-bonds are shown when only one of the two atoms is visible.","boolean","off","2" "hash_max","controls the maximum amount of memory used by the various temporary 3D hash tables used in geometry generation and raytracing.","integer","depends","0" "heavy_neighbor_cutoff","cutoff for finding heavy neighbors in measurements wizard.","float","3.5","0" diff --git a/layer1/SettingInfo.h b/layer1/SettingInfo.h index 5127b9004..f255d4aa8 100644 --- a/layer1/SettingInfo.h +++ b/layer1/SettingInfo.h @@ -902,6 +902,7 @@ enum { REC_f( 792, halogen_bond_as_acceptor_min_donor_angle , global , 120.0f ), REC_f( 793, halogen_bond_as_acceptor_min_acceptor_angle , global , 90.0f ), REC_f( 794, halogen_bond_as_acceptor_max_acceptor_angle , global , 170.0f ), + REC_f( 795, salt_bridge_distance , global , 5.0f ), #ifdef SETTINGINFO_IMPLEMENTATION #undef SETTINGINFO_IMPLEMENTATION diff --git a/layer2/ObjectDist.cpp b/layer2/ObjectDist.cpp index 2b9ff09cd..b7a7c7d1b 100644 --- a/layer2/ObjectDist.cpp +++ b/layer2/ObjectDist.cpp @@ -453,6 +453,9 @@ ObjectDist *ObjectDistNewFromSele(PyMOLGlobals * G, ObjectDist * oldObj, } else if (mode == 9) { // 9: halogen-bond interaction I->DSet[a].reset(pymol::FindHalogenBondInteractions(G, I->DSet[a].release(), sele1, state1, sele2, state2, cutoff, &dist)); + } else if (mode == 10) { // 10: salt-bridge interaction + I->DSet[a].reset(pymol::FindSaltBridgeInteractions(G, + I->DSet[a].release(), sele1, state1, sele2, state2, cutoff, &dist)); } else { I->DSet[a].reset(SelectorGetDistSet( G, I->DSet[a].release(), sele1, state1, sele2, state2, mode, cutoff, &dist)); diff --git a/layer3/Executive.cpp b/layer3/Executive.cpp index 8ec921ea0..a1795c25f 100644 --- a/layer3/Executive.cpp +++ b/layer3/Executive.cpp @@ -9539,7 +9539,10 @@ pymol::Result ExecutiveDistance(PyMOLGlobals* G, const char* nam, SettingSet(cSetting_dash_color, "0xff8800" /* light red */, obj); break; case 9: // "halogen-bonds" - SettingSet(cSetting_dash_color, "0xff00ff" /* magenta */, obj); + SettingSet(cSetting_dash_color, "0xaa00ff" /* dark magenta */, obj); + break; + case 10: // "salted-bridge" + SettingSet(cSetting_dash_color, "0xff2eff" /* light magenta */, obj); break; } diff --git a/layer3/Interactions.cpp b/layer3/Interactions.cpp index dde79cd26..d089a9daf 100644 --- a/layer3/Interactions.cpp +++ b/layer3/Interactions.cpp @@ -657,6 +657,12 @@ static bool CheckHalogenBondAsAcceptor(ObjectMolecule* don_obj, int don_atom, return result; } +SaltBridgeCriteria::SaltBridgeCriteria(PyMOLGlobals* G) +{ + m_distance = + SettingGet(G, nullptr, nullptr, cSetting_salt_bridge_distance); +} + HalogenBondCriteria::HalogenBondCriteria(PyMOLGlobals* G) { m_distance = @@ -673,25 +679,21 @@ HalogenBondCriteria::HalogenBondCriteria(PyMOLGlobals* G) cSetting_halogen_bond_as_acceptor_max_acceptor_angle); } -DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, - int state1, int sele2, int state2, float cutoff, float* result) +/** + * Prepare neighbor tables + * + * @param sele1 - selections index + * @param state1 - state index + * @param sele2 - selection index + * @param state2 - state index + * + * @return maximum number of atoms + */ +int PrepareNeighborTables( + PyMOLGlobals* G, int sele1, int state1, int sele2, int state2) { CSelector* I = G->Selector; - int exclusion = SettingGet(G, cSetting_h_bond_exclusion); - - int numVerts = 0; - *result = 0.0f; - // if the dist set exists, get info from it, otherwise get a new one - if (ds == nullptr) { - ds = DistSetNew(G); - } else { - numVerts = ds->NIndex; // number of vertices - } - - auto& coords = ds->Coord; - coords.reserve(10); - // update states: if the two are the same, update that one state, else update // all states if (state1 < 0 || state2 < 0 || state1 != state2) { @@ -727,6 +729,103 @@ DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, } } + return max_n_atom; +} + +/** + * Insert distance info into DistSet and copy atom's coordinates + * + * @param ds - distance set + * @param state1 - state index + * @param state2 - state index + * @param ai1 - atom info + * @param ai2 - atom info + * @param atom1_vv - atom coordinates + * @param atom2_vv - atom coordinates + * @param numVerts - number of vertices + */ +void InsertDistanceInfo(PyMOLGlobals* G, DistSet* ds, int state1, int state2, + AtomInfoType* ai1, AtomInfoType* ai2, float* atom1_vv, float* atom2_vv, + int numVerts) +{ + // Insert DistInfo records for updating distances + // Init/Add the elem to the DistInfo list + ds->MeasureInfo.emplace_front(); + auto* atom1Info = &ds->MeasureInfo.front(); + + // TH + atom1Info->id[0] = AtomInfoCheckUniqueID(G, ai1); + atom1Info->id[1] = AtomInfoCheckUniqueID(G, ai2); + + atom1Info->offset = numVerts; // offset into this DSet's Coord + atom1Info->state[0] = state1; // state1 of sel1 + atom1Info->state[1] = state2; + atom1Info->measureType = cRepDash; // DISTANCE-dash + + auto& coords = ds->Coord; + + // see if coords has room at another 6 floats + VLACheck(coords, float, (numVerts * 3) + 6); + float* vv0 = coords + (numVerts * 3); + + if (atom1_vv != nullptr && atom2_vv != nullptr) { + const size_t count_to_copy = 3; + + std::copy_n(atom1_vv, count_to_copy, vv0); + vv0 += count_to_copy; + atom1_vv += count_to_copy; + + std::copy_n(atom2_vv, count_to_copy, vv0); + vv0 += count_to_copy; + atom2_vv += count_to_copy; + } +} + +/** + * Create coverage vector + * + * @param sele1 - selection index + * @param sele2 - selection index + * @return vector of booleans which determines if a given atom appears in sele1 + * and sele2 + */ +static std::vector CreateCoverage(PyMOLGlobals* G, int sele1, int sele2) +{ + CSelector* I = G->Selector; + std::vector result(I->Table.size()); + + // coverage determines if a given atom appears in sele1 and sele2 + for (SelectorAtomIterator iter(I); iter.next();) { + int s = iter.getAtomInfo()->selEntry; + if (SelectorIsMember(G, s, sele1) && SelectorIsMember(G, s, sele2)) { + result[iter.a] = true; + } + } + + return result; +} + +DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, + int state1, int sele2, int state2, float cutoff, float* result) +{ + CSelector* I = G->Selector; + + int exclusion = SettingGet(G, cSetting_h_bond_exclusion); + + int numVerts = 0; + *result = 0.0f; + // if the dist set exists, get info from it, otherwise get a new one + if (ds == nullptr) { + ds = DistSetNew(G); + } else { + numVerts = ds->NIndex; // number of vertices + } + + auto& coords = ds->Coord; + coords.reserve(10); + + int max_n_atom = PrepareNeighborTables(G, sele1, state1, sele2, state2); + HalogenBondCriteria halogenbcRec(G); cutoff = halogenbcRec.m_distance; if (cutoff < 0.0f) { @@ -735,13 +834,7 @@ DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, } // coverage determines if a given atom appears in sel1 and sel2 - std::vector coverage(I->Table.size()); - for (SelectorAtomIterator iter(I); iter.next();) { - int s = iter.getAtomInfo()->selEntry; - if (SelectorIsMember(G, s, sele1) && SelectorIsMember(G, s, sele2)) { - coverage[iter.a] = true; - } - } + std::vector coverage = CreateCoverage(G, sele1, sele2); // this creates an interleaved list of ints for mapping ids to states within a // given neighborhood @@ -771,6 +864,10 @@ DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, int at1 = I->Table[a1].atom; int at2 = I->Table[a2].atom; + if (sele1 == sele2 && at1 > at2) { + continue; + } + // get the object for this global atom ID ObjectMolecule* obj1 = I->Obj[I->Table[a1].model]; ObjectMolecule* obj2 = I->Obj[I->Table[a2].model]; @@ -799,8 +896,6 @@ DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, // if we pass the bonding cutoff if (dist < cutoff) { - float h_crd[3]; - AtomInfoType* h_ai = nullptr; bool a_keeper = false; if (exclusion && obj1 == obj2) { @@ -846,44 +941,156 @@ DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, } } - if (sele1 == sele2 && at1 > at2) { - a_keeper = false; + if (a_keeper) { + InsertDistanceInfo( + G, ds, state1, state2, ai1, ai2, don_vv, acc_vv, numVerts); + dist_cnt++; + dist_sum += dist; + numVerts += 2; } + } + } + } + } + } + } - if (a_keeper) { - // Insert DistInfo records for updating distances - // Init/Add the elem to the DistInfo list - ds->MeasureInfo.emplace_front(); - auto* atom1Info = &ds->MeasureInfo.front(); + if (dist_cnt > 0) { + (*result) = dist_sum / dist_cnt; + } - // TH - atom1Info->id[0] = AtomInfoCheckUniqueID(G, ai1); - atom1Info->id[1] = AtomInfoCheckUniqueID(G, ai2); + if (coords) { + coords.resize((numVerts + 1) * 3); + } - atom1Info->offset = numVerts; // offset into this DSet's Coord - atom1Info->state[0] = state1; // state1 of sel1 - atom1Info->state[1] = state2; - atom1Info->measureType = cRepDash; // DISTANCE-dash + ds->NIndex = numVerts; - // we have a distance we want to keep - dist_cnt++; - dist_sum += dist; - // see if coords has room at another 6 floats - VLACheck(coords, float, (numVerts * 3) + 6); - float* vv0 = coords + (numVerts * 3); + return ds; +} + +DistSet* FindSaltBridgeInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, + int state1, int sele2, int state2, float cutoff, float* result) +{ + CSelector* I = G->Selector; + + int exclusion = SettingGet(G, cSetting_h_bond_exclusion); + + int numVerts = 0; + *result = 0.0f; + // if the dist set exists, get info from it, otherwise get a new one + if (ds == nullptr) { + ds = DistSetNew(G); + } else { + numVerts = ds->NIndex; // number of vertices + } - if (don_vv != nullptr && acc_vv != nullptr) { - const size_t count_to_copy = 3; + auto& coords = ds->Coord; + coords.reserve(10); - std::copy_n(don_vv, count_to_copy, vv0); - vv0 += count_to_copy; - don_vv += count_to_copy; + int max_n_atom = PrepareNeighborTables(G, sele1, state1, sele2, state2); - std::copy_n(acc_vv, count_to_copy, vv0); - vv0 += count_to_copy; - acc_vv += count_to_copy; + SaltBridgeCriteria saltbcRec(G); + cutoff = saltbcRec.m_distance; + if (cutoff < 0.0f) { + const float max_cutoff = 1000.0f; + cutoff = max_cutoff; + } + + // coverage determines if a given atom appears in sel1 and sel2 + std::vector coverage = CreateCoverage(G, sele1, sele2); + + // this creates an interleaved list of ints for mapping ids to states within a + // given neighborhood + std::vector interstate_vector = + SelectorGetInterstateVector(G, sele1, state1, sele2, state2, cutoff); + int cnt = interstate_vector.size() / 2; + + std::vector zero(max_n_atom); + std::vector scratch(max_n_atom); + + float dist_sum = 0.0f; + int dist_cnt = 0; + + // for each state + for (int a = 0; a < cnt; a++) { + + // get the interstate atom identifier for the two atoms to distance + int a1 = interstate_vector[a * 2]; + int a2 = interstate_vector[a * 2 + 1]; + + // check their coverage to avoid duplicates + if (a1 < a2 || (a1 != a2 && !(coverage[a1] && coverage[a2])) || + (state1 != state2)) { + // eliminate reverse duplicates + + // get the object-local atom ID + int at1 = I->Table[a1].atom; + int at2 = I->Table[a2].atom; + + if (sele1 == sele2 && at1 > at2) { + continue; + } + + // get the object for this global atom ID + ObjectMolecule* obj1 = I->Obj[I->Table[a1].model]; + ObjectMolecule* obj2 = I->Obj[I->Table[a2].model]; + + // the states are valid for these two atoms + if (state1 < obj1->NCSet && state2 < obj2->NCSet) { + // get the coordinate sets for both atoms + CoordSet* cs1 = obj1->CSet[state1]; + CoordSet* cs2 = obj2->CSet[state2]; + if (cs1 != nullptr && cs2 != nullptr) { + // for bonding + float* anion_vv = nullptr; + float* cation_vv = nullptr; + + // grab the appropriate atom information for this object-local atom + AtomInfoType* ai1 = obj1->AtomInfo + at1; + AtomInfoType* ai2 = obj2->AtomInfo + at2; + + int atom1_charge = ai1->formalCharge; + int atom2_charge = ai2->formalCharge; + if (atom1_charge * atom2_charge >= 0) { + continue; + } + + if (ai1->isHydrogen() || ai2->isHydrogen()) { + continue; + } + + int idx1 = cs1->atmToIdx(at1); + int idx2 = cs2->atmToIdx(at2); + + if (idx1 >= 0 && idx2 >= 0) { + // actual distance calculation from ptA to ptB + float dist = + (float) diff3f(cs1->coordPtr(idx1), cs2->coordPtr(idx2)); + + // if we pass the bonding cutoff + if (dist < cutoff) { + + bool a_keeper = false; + if (exclusion && obj1 == obj2) { + a_keeper = !SelectorCheckNeighbors( + G, exclusion, obj1, at1, at2, zero.data(), scratch.data()); + } + + if (a_keeper) { + if (atom1_charge < 0) { + anion_vv = cs1->coordPtr(idx1); + cation_vv = cs2->coordPtr(idx2); + } else { + anion_vv = cs2->coordPtr(idx2); + cation_vv = cs1->coordPtr(idx1); } + } + if (a_keeper) { + InsertDistanceInfo(G, ds, state1, state2, ai1, ai2, anion_vv, + cation_vv, numVerts); + dist_cnt++; + dist_sum += dist; numVerts += 2; } } diff --git a/layer3/Interactions.h b/layer3/Interactions.h index c92fd4cd4..505e250a9 100644 --- a/layer3/Interactions.h +++ b/layer3/Interactions.h @@ -46,6 +46,23 @@ class HalogenBondCriteria HalogenBondCriteria(PyMOLGlobals* G); }; +/** + * Salt Bridge criteria + */ +class SaltBridgeCriteria +{ + +public: + // Cutoff distance + float m_distance = 5.0f; + +public: + /** + * Construct a new Salt Bridge Criteria object and initialize it's member variables + */ + SaltBridgeCriteria(PyMOLGlobals* G); +}; + DistSet* FindPiInteractions(PyMOLGlobals* G, DistSet* ds, // int sele1, int state1, // @@ -68,4 +85,20 @@ DistSet* FindPiInteractions(PyMOLGlobals* G, */ DistSet* FindHalogenBondInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, int state1, int sele2, int state2, float cutoff, float* result); + +/** + * Find Salt-bridge interactions + * + * @param ds - DistSet + * @param sele1 - selections index + * @param state1 - state index + * @param sele2 - selection index + * @param state2 - state index + * @param cutoff - cutoff distance + * @param result - result average distance + * + * @return DistSet - distance set + */ +DistSet* FindSaltBridgeInteractions(PyMOLGlobals* G, DistSet* ds, int sele1, + int state1, int sele2, int state2, float cutoff, float* result); } diff --git a/modules/pymol/menu.py b/modules/pymol/menu.py index 101f4ed4a..eb59a806a 100644 --- a/modules/pymol/menu.py +++ b/modules/pymol/menu.py @@ -993,6 +993,14 @@ def halogen_bond(self_cmd, sele): ], ] +def salt_bridge(self_cmd, sele): + return [[2, 'Salt-Bridge Interactions:', '']] + [ + [ + 1, 'salt-bridge', 'cmd.distance("' + sele + '_salt_bridge","' + + sele + '","same",reset=1,mode=10)' + ], + ] + def polar(self_cmd, sele): return [[ 2, 'Polar Contacts:', ''], [ 1, 'within selection' , @@ -1047,6 +1055,7 @@ def find(self_cmd, sele): for d in (3.0, 3.5, 4.0) ]], [ 1, 'halogen-bond interactions', halogen_bond(self_cmd, sele)], + [ 1, 'salt-bridge interactions', salt_bridge(self_cmd, sele)], [ 1, 'pi interactions', [[ 2, 'Pi Interactions:', '']] + [ [1, 'all', 'cmd.pi_interactions("'+sele+'_pi_interactions","'+sele+'",reset=1)'], [1, 'pi-pi', 'cmd.distance("'+sele+'_pi_pi","'+sele+'","same",reset=1,mode=6)'],