Aromatic bond visualization

[dkoes]

Stick style should show aromatic bonds with solid/dashed representation.

[dxdc]

@dkoes

I've done more follow up on #639, and I have a working implementation that could also be used for this feature (have not pushed the PR yet). It's a more general solution in which one can specify partial bond order of any bond to control the thickness of the remaining dashed bond.

For example, instead of specifying a bond order of 2, it uses 1.7 in this example:

--

However, that brings up a few questions.

1. Is there a mathematical way to determine -- in a multi-bond situation -- which of the bonds should be dashed? Consider this example, with the dashed bonds reversed.

In my WIP implementation, I have done it in two ways, but a more mathematical default would be superior (e.g., to somehow prioritize the "shorter" bond as dashed).

• dashPriority property on the atom itself
• Bond order which ends with 001, e.g. 1.7001. This could be easier, in some cases, for tuning SDF files (particularly in the context of an animation) without trying to properly determine which atom it should be assigned to in the settings.

viewer.setStyle({serial: 4}, {
    stick: {
        radius: 0.2,
        dashPriority: true,
    }
});

2. In cases where atoms have different colors, what is the correct way to determine the dash color?

• Existing implementation (about half and half)

• Gradient color scheme, slowly transitioning from one to the other

• Calculate the most appropriate color by relative contrast, or another means

• Or, provide user option to override the color (perhaps, in addition)

[dxdc]

Here is another view of the coloration question (ignore the fact that this is a Cl in the example, it was just done for convenience). There may be situations that this kind of view is desirable, for example, an allylic transition state involving a carbonyl.

[dkoes]

I think the common convention is that the dash should be inside the ring. I personal prefer the discrete change in color halfway for the bond color.

@ghutchis Do you have any opinions as a molecular visualization expert?

[dxdc]

I think the common convention is that the dash should be inside the ring.

Of course. My question is: can that be determined in the context of available functions or some type of vector mathematics? It seems like the vector length between from/to are identical. It also seems like it could be a bit tricky to generally solve this for cases where the dashes are used to represent transition states.

EDIT: relevant section of code is here. We need to determine which cylinder should be drawn dashed.

3Dmol.js/src/GLModel.ts

Lines 906 to 955 in b038d7b

	if (atom.bondOrder[i] == 2) {
	r = bondR * doubleBondScale;
	v.multiplyScalar(r * 1.5);
	p1a = p1.clone();
	p1a.add(v);
	p1b = p1.clone();
	p1b.sub(v);
	p2a = p1a.clone();
	p2a.add(dir);
	p2b = p1b.clone();
	p2b.add(dir);
	if (C1 != C2) {
	mp = new Vector3().addVectors(p1a, p2a)
	.multiplyScalar(0.5);
	mp2 = new Vector3().addVectors(p1b, p2b)
	.multiplyScalar(0.5);
	drawCyl(geo, p1a, mp, r, C1, mfromCap, 0);
	drawCyl(geo, mp, p2a, r, C2, 0, mtoCap);
	drawCyl(geo, p1b, mp2, r, C1, mfromCap, 0);
	drawCyl(geo, mp2, p2b, r, C2, 0, mtoCap);
	} else {
	drawCyl(geo, p1a, p2a, r, C1, mfromCap, mtoCap);
	drawCyl(geo, p1b, p2b, r, C1, mfromCap, mtoCap);
	}
	atomneedsi = atom.clickable || atom.hoverable;
	atom2needsi = atom2.clickable || atom2.hoverable;
	if (atomneedsi || atom2needsi) {
	if (!mp) mp = new Vector3().addVectors(p1a, p2a)
	.multiplyScalar(0.5);
	if (!mp2) mp2 = new Vector3().addVectors(p1b, p2b)
	.multiplyScalar(0.5);
	if (atomneedsi) {
	cylinder1a = new Cylinder(p1a, mp, r);
	cylinder1b = new Cylinder(p1b, mp2, r);
	atom.intersectionShape.cylinder.push(cylinder1a);
	atom.intersectionShape.cylinder.push(cylinder1b);
	}
	if (atom2needsi) {
	cylinder2a = new Cylinder(p2a, mp, r);
	cylinder2b = new Cylinder(p2b, mp2, r);
	atom2.intersectionShape.cylinder.push(cylinder2a);
	atom2.intersectionShape.cylinder.push(cylinder2b);
	}
	}
	}

[dkoes]

I think you should be able to get it from v which is a function of the larger atomic environment. If it doesn't consistently point in the concave direction, we can probably make it (but I don't have time to look at it in detail now).

[dxdc]

get it from v

Does that mean the order here in the code (cylinder 1: p1a/p2a and cylinder 2: p1b/p2b) should be consistent, or that there is some other calculation needed?

The existing order is not entirely consistent. This is an example from a Diels-Alder transition state.

Ultimately, in my implementation, I defined the logic like this. If there is a way to determine isSolidBondFirst in a consistent way, it could obviate the need for custom params to that end.

if (isSolidBondFirst || r === r2) {
	drawCylSolid(geo, p1a, p2a, ...
	drawCyl(geo, p1b, p2b, ...
} else {
	drawCyl(geo, p1a, p2a, ...
	drawCylSolid(geo, p1b, p2b, ...
}

[dxdc]

@dkoes, any update on this? Appreciate your help!

[dkoes]

Sorry for delay - life and work.

The most relevant function is getGoodCross - this is called from getSideBondV to pick the vector to define the plane the double bond is placed in. The current criteria is to find the "most divergent neighbor". This is probably resulting in a hydrogen being picked in the above example. I think if you changed to logic to prefer non-hydrogens (as long as the chosen atom isn't colinear) that would fix the transition state picture. Another heuristic is to prefer atoms connected with an aromatic bound. The idea is for v to always point in a consistent direction when part of a ring system.

Can you provide the structure for the transition state above? Definitely want to see this added as a test case when the WIP is done.

[dxdc]

thanks for the reply @dkoes. Maybe I don't understand how to use v.

consider the case where the colors are identical, e.g.

                            drawCyl(geo, p1a, p2a, r, C1, mfromCap, mtoCap);
                            drawCyl(geo, p1b, p2b, r, C1, mfromCap, mtoCap);

how do I decide (from v) which bond to draw dashed?

                    if (atom.bondOrder[i] == 2) {
                        r = bondR * doubleBondScale;


                        v.multiplyScalar(r * 1.5);
                        p1a = p1.clone();
                        p1a.add(v);
                        p1b = p1.clone();
                        p1b.sub(v);


                        p2a = p1a.clone();
                        p2a.add(dir);
                        p2b = p1b.clone();
                        p2b.add(dir);


                        if (C1 != C2) {
                            mp = new Vector3().addVectors(p1a, p2a)
                                .multiplyScalar(0.5);
                            mp2 = new Vector3().addVectors(p1b, p2b)
                                .multiplyScalar(0.5);
                            drawCyl(geo, p1a, mp, r, C1, mfromCap, 0);
                            drawCyl(geo, mp, p2a, r, C2, 0, mtoCap);
                            drawCyl(geo, p1b, mp2, r, C1, mfromCap, 0);
                            drawCyl(geo, mp2, p2b, r, C2, 0, mtoCap);
                        } else {
                            drawCyl(geo, p1a, p2a, r, C1, mfromCap, mtoCap);
                            drawCyl(geo, p1b, p2b, r, C1, mfromCap, mtoCap);
                        }

[dkoes]

If v is consistent, then p1a will either always be inside the ring or outside - so you should always draw p1a-p2a dashed and adjust the setting of v so that is the right answer (or always draw p1b-p2b).

[dxdc]

always draw p1b-p2b

thanks @dkoes that was very helpful! using that insight I made some more progress. benzene is drawn properly now, and my transition state is a lot closer (2/3 bonds are correct). note that this is without adjusting the getGoodCross function. I had initially tried to alter it, but had mixed results.

I tried to alter getGoodCross to compute a score, something like:

var score = (atom3.elem !== 'H' ? 1 : 0) + atom.bondOrder[j];

but I'm not sure how it may impact the existing calculations, and if the length should still be trusted or used. One interesting thing is that if I comment this line (e.g., to no longer return early), it actually returns a very strange result in that case -- essentially updating the bond orientation out of plane by 90 degrees. It looks pretty close with the current getGoodCross implementation, so I'm hesitant to add too many other conditions, but welcome any feedback you have!

                        if(l > bestlen) {
                            bestlen = l;
                            bestv = v;
                            if(bestlen > 0.1) {
                                // return bestv;  // comment this line to no longer return early
                            }
                        }

Now, double-bond is rotated 90 degrees for some reason:

SDF file:

M0001

Created by GaussView 6.0.16
 16 16  0  0  0  0  0  0  0  0  0    0
   -1.5030   -0.6101   -0.7999 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.5784   -0.5027   -0.7480 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.0924   -1.4430   -0.2494 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7456    0.2411   -1.4556 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2085    1.0672   -1.9846 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.7456    0.2411   -1.4556 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.2085    1.0672   -1.9846 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.5030   -0.6101   -0.7999 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.0924   -1.4430   -0.2494 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.5784   -0.5027   -0.7480 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6530    0.4190    1.7575 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.2284    1.2692    1.4141 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.2299   -0.4405    2.0656 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6530    0.4190    1.7575 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2284    1.2692    1.4141 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2299   -0.4405    2.0656 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1.5  0  0  0  0
  1 14  1  0  0  0  0
  4  5  1  0  0  0  0
  4  6  4  0  0  0  0
  6  7  1  0  0  0  0
  6  8  1.5  0  0  0  0
  8  9  1  0  0  0  0
  8 10  1  0  0  0  0
  8 11  1  0  0  0  0
 11 12  1  0  0  0  0
 11 13  1  0  0  0  0
 11 14  1.5  0  0  0  0
 14 15  1  0  0  0  0
 14 16  1  0  0  0  0
M  END
$$$$

[dkoes]

The double bond will be in the plane defined by atom, atom2 and the atom3 that produces the bestv. You are getting a rotation because the hydrogen is being selected as the best atom3. I would skip the score and simply ignore hydrogens unless no non-hydrogen atoms produce a large enough bestv.

[dkoes]

So the logic for picking the best atom might something like:
Is the previous best length < 0.1 and the current > 0.1? update best
Is the previous best a hydrogen and the current not (and length > 0.1)? update best
Is the previous best not an aromatic bound from atom2 to atom3 and the current is (and length > 0.1)? update best
Is the previous best the same in terms of hydrogen/aromatic but the length is larger? update best

[dxdc]

thanks for the continued discussion @dkoes.

As opposed to iterating through the atoms in getGoodCross and potentially returning early, I just indexed each one like this. My thought was that we could sort the candidates array afterwards based on a series of criteria.

I also included the total number of bonds in atom3 -- thought that could be interesting to determine a 'divergent' candidate. I would prefer to use something more quantitative than just "aromatic", if that is possible, since a) I'm not sure how to determine aromatic bound from atom2 to atom3 and b) I'm modifying bondOrder anyway in the SDF file to achieve a certain look.

candidates.push({
	v: v,
	l: l,
	nonHydrogen: atom3.elem !== 'H',
	numBonds: atom3.bonds.length
});

// ...then, later:

// Sort candidates: non-hydrogen > numBonds > longer length
candidates.sort((a, b) => {
	if (a.isNonHydrogen !== b.isNonHydrogen) {
		return a.isNonHydrogen ? -1 : 1;
	}

	if (a.numBonds !== b.numBonds) {
		 return b.numBonds - a.numBonds;
	}
				   
	return b.l - a.l;
});

After computing these properties, I outputted the candidates array and determined/selected the correct choice empirically, and have indicated it as such below (this is an unsorted list).

I'm wondering if part of the problem is the (perhaps) arbitrary selection of atom vs. atom2? If I reverse the choice of atom vs. atom2 in this case, it seems to fix it. But, this doesn't appear to work generally...e.g.., it yields an unexpected result for a benzene ring with a carbonyl.

if (atom2.bonds.length < atom.bonds.length) {
    [atom, atom2] = [atom2, atom];
}

The first entry there is giving a result that doesn't seem to match with your specs. For reference, that is the bond in the top left of this image. If it's helpful, from the SDF file I posted above, the relevant atoms (i.e., [atom.serial, atom2.serial, atom3.serial]) are:

• 0: [0, 3, 1]
• 1: [0, 3, 2] **this one yields the correctly rendered result
• 2: [0, 3, 13]

[dkoes]

Do you think you have a better goodCross?

I feel like the overall logic of choosing the bond plane and direction should be simplifiable, but I keep getting stuck with various corner cases when I try.

[dxdc]

No, I had the same experience and am not familiar enough with the codebase to proceed any further. I'm wondering if the logic for this is already solved by some open source Chemistry drawing app.

[dkoes]

Let's not make the perfect the enemy of the good - the aromatic dashed bonds would still be great to have. Can you make a PR? We can make it optional.

[dkoes]

Have you verified that your weird molecules above visualize the way you expect in other software? Maybe you are setting a standard higher than what the community currently accepts. If they look right in an open source viewer I can dig into the source code (e.g. JMol or pymol).

[dxdc]

@dkoes sorry for the late reply! It has been all hands on deck dealing with various stages of the plague here 😅

Have you verified that your weird molecules above visualize the way you expect in other software?

Hmm, you may be right. It's wrong there too, but it does look more consistent in pymol though (didn't test jmol). At least two of the bonds are out of plane, but interestingly, they are consistently off whereas I'm getting more of a mixed behavior here. I am still wondering if the distinction between atom vs. atom2 needs to be picked in some kind of more sophisticated way to get the desired behavior.

Unfortunately, my code to trial all this is unpolished and hacky, utilizing a specialty file format and not really suitable for distribution/PR unfortunately. I'll see if I can find time to keep working on this!

[dkoes]

My general feeling is that a solution that works for normal aromatic rings (flat, bonds to hydrogens are shorter than the aromatic bonds) is good enough.