in layer2/ObjectMolecule2.cpp [4064:4409]
int ObjectMoleculeConnect(ObjectMolecule * I, int *nbond, BondType ** bond, AtomInfoType * ai,
struct CoordSet *cs, int bondSearchMode,
int connectModeOverride)
{
#define cMULT 1
PyMOLGlobals *G = I->Obj.G;
int a, b, c, d, e, f, i, j;
int a1, a2;
float *v1, *v2;
int maxBond;
MapType *map;
int nBond;
BondType *ii1, *ii2;
int flag;
int order;
AtomInfoType *ai1, *ai2;
/* Sulfur cutoff */
float cutoff_s;
float cutoff_v;
float max_cutoff;
int repeat = true;
int discrete_chains = SettingGetGlobal_i(G, cSetting_pdb_discrete_chains);
int connect_bonded = SettingGetGlobal_b(G, cSetting_connect_bonded);
int connect_mode = SettingGetGlobal_i(G, cSetting_connect_mode);
int unbond_cations = SettingGetGlobal_i(G, cSetting_pdb_unbond_cations);
int ok = true;
cutoff_v = SettingGetGlobal_f(G, cSetting_connect_cutoff);
cutoff_s = cutoff_v + 0.2F;
max_cutoff = cutoff_s;
if(connectModeOverride >= 0)
connect_mode = connectModeOverride;
if(connect_mode == 2) { /* force use of distance-based connectivity,
ignoring that provided with file */
bondSearchMode = true;
cs->NTmpBond = 0;
VLAFreeP(cs->TmpBond);
} else if (connect_mode == 4) {
// mmCIF specific, fall back to default to get any bonds for PDB, XYZ, etc.
connect_mode = 0;
}
/* FeedbackMask[FB_ObjectMolecule]=0xFF; */
nBond = 0;
maxBond = cs->NIndex * 8;
(*bond) = VLACalloc(BondType, maxBond);
CHECKOK(ok, (*bond));
while(ok && repeat) {
repeat = false;
/*
* BOND SEARCH MODE
*/
if(cs->NIndex && bondSearchMode) { /* &&(!I->DiscreteFlag) WLD 010527 */
/* if there are atoms, and we need to search for bonds, instead of using
* (possibly) supplied CONECT records... */
PRINTFB(G, FB_ObjectMolecule, FB_Blather)
" ObjectMoleculeConnect: Searching for bonds amongst %d coordinates.\n",
cs->NIndex ENDFB(G);
if(Feedback(G, FB_ObjectMolecule, FB_Debugging)) {
for(a = 0; a < cs->NIndex; a++)
printf(" ObjectMoleculeConnect: coord %d %8.3f %8.3f %8.3f\n",
a, cs->Coord[a * 3], cs->Coord[a * 3 + 1], cs->Coord[a * 3 + 2]);
}
switch (connect_mode) {
/* DISTANCE BASED CONNECTIONS */
case 0: /* distance-based and explicit (not HETATM to HETATM) */
case 3: /* distance-based and explicit (even HETATM to HETATM) */
case 2: /* distance-based only */ {
/* distance-based bond location */
int violations = 0;
int *cnt = Alloc(int, cs->NIndex);
int valcnt;
CHECKOK(ok, cnt);
if (ok){
/* initialize each atom's (max) expected valence */
for(i = 0; i < cs->NIndex; i++) {
valcnt = AtomInfoGetExpectedValence(G, ai + cs->IdxToAtm[i]);
if(valcnt < 0)
valcnt = 6;
cnt[i] = valcnt;
}
}
/* make a map of the local neighborhood in space */
if (ok)
map = MapNew(G, max_cutoff + MAX_VDW, cs->Coord, cs->NIndex, NULL);
CHECKOK(ok, map);
if(ok) {
int dim12 = map->D1D2;
int dim2 = map->Dim[2];
for(i = 0; ok && i < cs->NIndex; i++) {
if(nBond > maxBond)
break;
/* atom i's position in space */
v1 = cs->Coord + (3 * i);
a1 = cs->IdxToAtm[i];
ai1 = ai + a1;
MapLocus(map, v1, &a, &b, &c);
/* d = [a-1, a, a+1] */
for(d = a - 1; ok && d <= a + 1; d++) {
int *j_ptr1 = map->Head + d * dim12 + (b - 1) * dim2;
/* e = [b-1, b, b+1] */
for(e = b - 1; ok && e <= b + 1; e++) {
int *j_ptr2 = j_ptr1 + c - 1;
j_ptr1 += dim2;
/* f = [c-1, c, c+1] */
for(f = c - 1; ok && f <= c + 1; f++) {
j = *(j_ptr2++); /* *MapFirst(map,d,e,f)); */
while(ok && j >= 0) {
if(i < j) {
/* position in space for atom 2 */
v2 = cs->Coord + (3 * j);
a2 = cs->IdxToAtm[j];
ai2 = ai + a2;
flag = is_distance_bonded(G, cs, ai1, ai2, v1, v2,
cutoff_v, connect_mode, discrete_chains,
connect_bonded, unbond_cations);
{
/* we have a bond, now process it */
if(flag) {
VLACheck((*bond), BondType, nBond);
CHECKOK(ok, (*bond));
if (ok){
(*bond)[nBond].index[0] = a1;
(*bond)[nBond].index[1] = a2;
(*bond)[nBond].stereo = 0;
order = 1;
}
/* if we allow bonds between chains and it screws up the
* bonding, disallow inter-chain bonds */
if(ok && discrete_chains < 0) { /* if we're allowing bonds between chains,
then make sure things don't get out of hand */
if(cnt[i] == -1)
violations++;
if(cnt[j] == -1)
violations++;
/* decrement free valences, since we have a bond */
cnt[i]--;
cnt[j]--;
if(violations > (cs->NIndex >> 3)) {
/* if more than 12% of the structure has excessive #'s of bonds... */
PRINTFB(G, FB_ObjectMolecule, FB_Blather)
" ObjectMoleculeConnect: Assuming chains are discrete...\n"
ENDFB(G);
discrete_chains = 1;
repeat = true;
goto do_it_again;
}
}
if(!ai1->hetatm || ai1->resn == G->lex_const.MSE) {
if(AtomInfoSameResidue(I->Obj.G, ai1, ai2)) {
/* hookup standard disconnected PDB residue */
assign_pdb_known_residue(G, ai1, ai2, &order);
}
}
(*bond)[nBond].order = -order; /* store tentative valence as negative */
nBond++;
}
}
}
j = MapNext(map, j);
}
}
}
}
}
do_it_again:
MapFree(map);
FreeP(cnt);
}
}
case 1: /* only use explicit connectivity from file (don't do anything) */
break;
case 4: /* dictionary-based connectivity */
/* TODO */
break;
}
PRINTFB(G, FB_ObjectMolecule, FB_Blather)
" ObjectMoleculeConnect: Found %d bonds.\n", nBond ENDFB(G);
if(Feedback(G, FB_ObjectMolecule, FB_Debugging)) {
for(a = 0; a < nBond; a++)
printf(" ObjectMoleculeConnect: bond %d ind0 %d ind1 %d\n",
a, (*bond)[a].index[0], (*bond)[a].index[1]);
}
}
if(repeat) {
nBond = 0;
}
}
/* if we have explicit connectivity, determine if we need to set check_conect_all */
if(ok && cs->NTmpBond && cs->TmpBond) {
int check_conect_all = false;
int pdb_conect_all = false;
PRINTFB(G, FB_ObjectMolecule, FB_Blather)
" ObjectMoleculeConnect: incorporating explicit bonds. %d %d\n",
nBond, cs->NTmpBond ENDFB(G);
if((nBond == 0) && (cs->NTmpBond > 0) &&
bondSearchMode && (connect_mode == 0) && cs->NIndex) {
/* if no bonds were found, and we have explicit connectivity,
* try to determine if we need to set pdb_conect_mode */
for(i = 0; i < cs->NIndex; i++) {
a1 = cs->IdxToAtm[i];
ai1 = ai + a1;
if(ai1->bonded && (!ai1->hetatm)) {
/* apparent PDB ATOM record with explicit bonding... */
check_conect_all = true;
break;
}
}
}
VLACheck((*bond), BondType, (nBond + cs->NTmpBond));
CHECKOK(ok, (*bond));
if (ok){
ii1 = (*bond) + nBond;
ii2 = cs->TmpBond;
}
if (ok) {
int n_atom = I->NAtom;
for(a = 0; a < cs->NTmpBond; a++) {
a1 = cs->IdxToAtm[ii2->index[0]]; /* convert bonds from index space */
a2 = cs->IdxToAtm[ii2->index[1]]; /* to atom space */
if((a1 >= 0) && (a2 >= 0) && (a1 < n_atom) && (a2 < n_atom)) {
if(check_conect_all) {
if((!ai[a1].hetatm) && (!ai[a2].hetatm)) {
/* found bond between non-HETATMs -- so tell PyMOL to CONECT all ATOMs
* when writing out a PDB file */
pdb_conect_all = true;
}
}
ai[a1].bonded = true;
ai[a2].bonded = true;
(*ii1) = (*ii2); /* note this copies owned ids and thus settings etc. */
ii1->index[0] = a1;
ii1->index[1] = a2;
ii1++;
ii2++;
}
}
}
if (ok){
nBond = nBond + cs->NTmpBond;
VLAFreeP(cs->TmpBond);
cs->NTmpBond = 0;
if(pdb_conect_all) {
int dummy;
if(!SettingGetIfDefined_b(G, I->Obj.Setting, cSetting_pdb_conect_all, &dummy)) {
CSetting **handle = NULL;
if(I->Obj.fGetSettingHandle) {
handle = I->Obj.fGetSettingHandle(&I->Obj, -1);
if(handle) {
SettingCheckHandle(G, handle);
SettingSet_b(*handle, cSetting_pdb_conect_all, true);
}
}
}
}
}
}
/* Link b/t ligand and residue? */
if(ok && cs->NTmpLinkBond && cs->TmpLinkBond) {
PRINTFB(G, FB_ObjectMolecule, FB_Blather)
"ObjectMoleculeConnect: incorporating linkage bonds. %d %d\n",
nBond, cs->NTmpLinkBond ENDFB(G);
VLACheck((*bond), BondType, (nBond + cs->NTmpLinkBond));
CHECKOK(ok, (*bond));
if (ok){
ii1 = (*bond) + nBond;
ii2 = cs->TmpLinkBond;
for(a = 0; a < cs->NTmpLinkBond; a++) {
a1 = ii2->index[0]; /* first atom is in object */
a2 = cs->IdxToAtm[ii2->index[1]]; /* second is in the cset */
ai[a1].bonded = true;
ai[a2].bonded = true;
(*ii1) = (*ii2); /* note this copies owned ids and thus settings etc. */
ii1->index[0] = a1;
ii1->index[1] = a2;
ii1++;
ii2++;
}
nBond = nBond + cs->NTmpLinkBond;
}
VLAFreeP(cs->TmpLinkBond);
cs->NTmpLinkBond = 0;
}
PRINTFD(G, FB_ObjectMolecule)
" ObjectMoleculeConnect: elminating duplicates with %d bonds...\n", nBond ENDFD;
if(ok && !I->DiscreteFlag) {
UtilSortInPlace(G, (*bond), nBond, sizeof(BondType), (UtilOrderFn *) BondInOrder);
if(nBond) { /* eliminate duplicates */
ii1 = (*bond) + 1;
ii2 = (*bond) + 1;
a = nBond - 1;
nBond = 1;
if(a > 0)
while(a--) {
if((ii2->index[0] != (ii1 - 1)->index[0]) ||
(ii2->index[1] != (ii1 - 1)->index[1])) {
*(ii1++) = *(ii2++); /* copy bond */
nBond++;
} else {
if((ii2->order > 0) && ((ii1 - 1)->order < 0))
(ii1 - 1)->order = ii2->order; /* use most certain valence */
ii2++; /* skip bond */
}
}
VLASize((*bond), BondType, nBond);
CHECKOK(ok, (*bond));
}
}
/* restore bond oder positivity */
if (ok){
ii1 = *bond;
for(a = 0; a < nBond; a++) {
if(ii1->order < 0)
ii1->order = -ii1->order;
ii1++;
}
}
PRINTFD(G, FB_ObjectMolecule)
" ObjectMoleculeConnect: leaving with %d bonds...\n", nBond ENDFD;
if (ok)
*nbond = nBond;
return ok;
}