wip, bond order rendering

src/mol-geo/representation/structure/ball-and-stick.ts

@@ -18,7 +18,7 @@ export const DefaultBallAndStickProps = {
     ...DefaultElementSphereProps,
     ...DefaultIntraUnitLinkProps,
 
-    sizeTheme: { name: 'physical', factor: 0.2 } as SizeTheme,
+    sizeTheme: { name: 'uniform', value: 0.25 } as SizeTheme,
 }
 export type BallAndStickProps = Partial<typeof DefaultBallAndStickProps>
 
@@ -31,7 +31,7 @@ export function BallAndStickRepresentation(): StructureRepresentation<BallAndSti
             return [ ...sphereRepr.renderObjects, ...intraLinkRepr.renderObjects ]
         },
         create: (structure: Structure, props: BallAndStickProps = {} as BallAndStickProps) => {
-            const p = Object.assign({}, props, DefaultBallAndStickProps)
+            const p = Object.assign({}, DefaultBallAndStickProps, props)
             return Task.create('Creating BallAndStickRepresentation', async ctx => {
                 await sphereRepr.create(structure, p).runInContext(ctx)
                 await intraLinkRepr.create(structure, p).runInContext(ctx)
@@ -49,7 +49,7 @@ export function BallAndStickRepresentation(): StructureRepresentation<BallAndSti
             if (isEmptyLoci(sphereLoci)) {
                 return intraLinkLoci
             } else {
-                return sphereLoci 
+                return sphereLoci
             }
         },
         mark: (loci: Loci, action: MarkerAction) => {

src/mol-geo/representation/structure/visual/intra-unit-link-cylinder.ts

@@ -28,42 +28,121 @@ import { MarkerAction, applyMarkerAction, createMarkers, MarkerData } from '../.
 import { SizeTheme } from '../../../theme';
 import { chainIdLinkColorData } from '../../../theme/structure/color/chain-id';
 
-async function createLinkCylinderMesh(ctx: RuntimeContext, unit: Unit, mesh?: Mesh) {
+const DefaultLinkCylinderProps = {
+    linkScale: 0.4,
+    linkSpacing: 1,
+    linkRadius: 0.25,
+    radialSegments: 16
+}
+type LinkCylinderProps = typeof DefaultLinkCylinderProps
+
+async function createLinkCylinderMesh(ctx: RuntimeContext, unit: Unit, props: LinkCylinderProps, mesh?: Mesh) {
     if (!Unit.isAtomic(unit)) return Mesh.createEmpty(mesh)
 
     const elements = unit.elements;
     const links = unit.links
-    const { edgeCount, a, b } = links
+    const { edgeCount, a, b, edgeProps, offset } = links
+    const orders = edgeProps.order
 
     if (!edgeCount) return Mesh.createEmpty(mesh)
 
-    // TODO calculate vertextCount properly
+    // approximate vertextCount, exact calculation would need to take link orders into account
     const vertexCount = 32 * edgeCount
     const meshBuilder = MeshBuilder.create(vertexCount, vertexCount / 2, mesh)
 
     const va = Vec3.zero()
     const vb = Vec3.zero()
-    const vt = Vec3.zero()
+    const vd = Vec3.zero()
+    const vc = Vec3.zero()
     const m = Mat4.identity()
+    const mt = Mat4.identity()
+
+    const vShift = Vec3.zero()
+    const vCenter = Vec3.zero()
+    const vRef = Vec3.zero()
+
+    const { linkScale, linkSpacing, linkRadius, radialSegments } = props
+
+    const cylinderParams = {
+        height: 1,
+        radiusTop: linkRadius,
+        radiusBottom: linkRadius,
+        radialSegments,
+        openEnded: true
+    }
 
     const pos = unit.conformation.invariantPosition
     // const l = Element.Location()
     // l.unit = unit
 
+    // assumes aI < bI
+    function getRefPos(aI: number, bI: number) {
+        if (aI > bI) console.log('aI > bI')
+        for (let i = offset[aI], il = offset[aI + 1]; i < il; ++i) {
+            if (b[i] !== bI) return pos(elements[b[i]], vRef)
+        }
+        for (let i = offset[bI], il = offset[bI + 1]; i < il; ++i) {
+            if (a[i] !== aI) return pos(elements[a[i]], vRef)
+        }
+        // console.log('no ref', aI, bI, unit.model.atomicHierarchy.atoms.auth_atom_id.value(aI), unit.model.atomicHierarchy.atoms.auth_atom_id.value(bI), offset[aI], offset[aI + 1], offset[bI], offset[bI + 1], offset)
+        return null
+    }
+
     for (let edgeIndex = 0, _eI = edgeCount * 2; edgeIndex < _eI; ++edgeIndex) {
         const aI = elements[a[edgeIndex]], bI = elements[b[edgeIndex]];
-        // each edge is included twice because of the "adjacency list" structure
-        // keep only the 1st occurence.
-        if (aI >= bI) continue;
+        // Each edge is included twice to allow for coloring/picking
+        // the half closer to the first vertex, i.e. vertex a.
         pos(aI, va)
         pos(bI, vb)
+        const d = Vec3.distance(va, vb)
+
+        Vec3.sub(vd, vb, va)
+        Vec3.scale(vd, Vec3.normalize(vd, vd), d / 4)
+        Vec3.add(vc, va, vd)
+        // ensure both edge halfs are pointing in the the same direction so the triangles align
+        if (aI > bI) Vec3.scale(vd, vd, -1)
+        Vec3.makeRotation(m, Vec3.create(0, 1, 0), vd)
 
-        Vec3.scale(vt, Vec3.add(vt, va, vb), 0.5)
-        Vec3.makeRotation(m, Vec3.create(0, 1, 0), Vec3.sub(vb, vb, va))
-        Mat4.setTranslation(m, vt)
-        
+        const order = orders[edgeIndex]
         meshBuilder.setId(edgeIndex)
-        meshBuilder.addCylinder(m, { radiusTop: 0.2, radiusBottom: 0.2 })
+        cylinderParams.height = d / 2
+
+        if (order === 2 || order === 3) {
+            const multiRadius = linkRadius * (linkScale / (0.5 * order))
+            const absOffset = (linkRadius - multiRadius) * linkSpacing
+
+            if (aI < bI) {
+                calculateShiftDir(vShift, va, vb, getRefPos(a[edgeIndex], b[edgeIndex]))
+            } else {
+                calculateShiftDir(vShift, vb, va, getRefPos(b[edgeIndex], a[edgeIndex]))
+            }
+            Vec3.setMagnitude(vShift, vShift, absOffset)
+
+            cylinderParams.radiusTop = multiRadius
+            cylinderParams.radiusBottom = multiRadius
+
+            if (order === 3) {
+                Mat4.fromTranslation(mt, vc)
+                Mat4.mul(mt, mt, m)
+                meshBuilder.addCylinder(mt, cylinderParams)
+            }
+
+            Vec3.add(vCenter, vc, vShift)
+            Mat4.fromTranslation(mt, vCenter)
+            Mat4.mul(mt, mt, m)
+            meshBuilder.addCylinder(mt, cylinderParams)
+
+            Vec3.sub(vCenter, vc, vShift)
+            Mat4.fromTranslation(mt, vCenter)
+            Mat4.mul(mt, mt, m)
+            meshBuilder.addCylinder(mt, cylinderParams)
+        } else {
+            cylinderParams.radiusTop = linkRadius
+            cylinderParams.radiusBottom = linkRadius
+
+            Mat4.setTranslation(m, vc)
+            meshBuilder.addCylinder(m, cylinderParams)
+        }
 
         if (edgeIndex % 10000 === 0 && ctx.shouldUpdate) {
             await ctx.update({ message: 'Cylinder mesh', current: edgeIndex, max: edgeCount });
@@ -75,6 +154,7 @@ async function createLinkCylinderMesh(ctx: RuntimeContext, unit: Unit, mesh?: Me
 
 export const DefaultIntraUnitLinkProps = {
     ...DefaultStructureProps,
+    ...DefaultLinkCylinderProps,
     sizeTheme: { name: 'physical', factor: 0.3 } as SizeTheme,
     flipSided: false,
     flatShaded: false,
@@ -87,7 +167,6 @@ export function IntraUnitLinkVisual(): UnitsVisual<IntraUnitLinkProps> {
     let currentProps: typeof DefaultIntraUnitLinkProps
     let mesh: Mesh
     let currentGroup: Unit.SymmetryGroup
-    // let vertexMap: VertexMap
 
     return {
         renderObjects,
@@ -101,7 +180,7 @@ export function IntraUnitLinkVisual(): UnitsVisual<IntraUnitLinkProps> {
             const elementCount = Unit.isAtomic(unit) ? unit.links.edgeCount * 2 : 0
             const instanceCount = group.units.length
 
-            mesh = await createLinkCylinderMesh(ctx, unit)
+            mesh = await createLinkCylinderMesh(ctx, unit, currentProps)
 
             if (ctx.shouldUpdate) await ctx.update('Computing link transforms');
             const transforms = createTransforms(group)
@@ -201,7 +280,7 @@ function markLink(loci: Loci, action: MarkerAction, group: Unit.SymmetryGroup, v
         for (const b of loci.links) {
             const unitIdx = Unit.findUnitById(b.aUnit.id, group.units)
             if (unitIdx !== -1) {
-                const _idx = unit.links.getEdgeIndex(b.aIndex, b.bIndex)
+                const _idx = unit.links.getDirectedEdgeIndex(b.aIndex, b.bIndex)
                 if (_idx !== -1) {
                     const idx = _idx
                     if (applyMarkerAction(array, idx, idx + 1, action) && !changed) {
@@ -216,4 +295,34 @@ function markLink(loci: Loci, action: MarkerAction, group: Unit.SymmetryGroup, v
     if (changed) {
         ValueCell.update(tMarker, tMarker.ref.value)
     }
+}
+
+const tmpShiftV12 = Vec3.zero()
+const tmpShiftV13 = Vec3.zero()
+
+/** Calculate 'shift' direction that is perpendiculat to v1 - v2 and goes through v3 */
+function calculateShiftDir (out: Vec3, v1: Vec3, v2: Vec3, v3: Vec3 | null) {
+    Vec3.sub(tmpShiftV12, v1, v2)
+
+    if (v3 !== null) {
+        Vec3.sub(tmpShiftV13, v1, v3)
+    } else {
+        Vec3.copy(tmpShiftV13, v1)  // no reference point, use v1
+    }
+    Vec3.normalize(tmpShiftV13, tmpShiftV13)
+
+    // ensure v13 and v12 are not colinear
+    let dp = Vec3.dot(tmpShiftV12, tmpShiftV13)
+    if (1 - Math.abs(dp) < 1e-5) {
+        Vec3.set(tmpShiftV13, 1, 0, 0)
+        dp = Vec3.dot(tmpShiftV12, tmpShiftV13)
+        if (1 - Math.abs(dp) < 1e-5) {
+            Vec3.set(tmpShiftV13, 0, 1, 0)
+            dp = Vec3.dot(tmpShiftV12, tmpShiftV13)
+        }
+    }
+
+    Vec3.setMagnitude(tmpShiftV12, tmpShiftV12, dp)
+    Vec3.sub(tmpShiftV13, tmpShiftV13, tmpShiftV12)
+    return Vec3.normalize(out, tmpShiftV13)
 }
\ No newline at end of file

src/mol-math/graph/int-adjacency-graph.ts

@@ -28,8 +28,17 @@ interface IntAdjacencyGraph<EdgeProps extends IntAdjacencyGraph.EdgePropsBase =
      *
      * Because the a and b arrays contains each edge twice,
      * this always returns the smaller of the indices.
+     *
+     * `getEdgeIndex(i, j) === getEdgeIndex(j, i)`
      */
     getEdgeIndex(i: number, j: number): number,
+    /**
+     * Get the edge index between i-th and j-th vertex.
+     * -1 if the edge does not exist.
+     *
+     * `getEdgeIndex(i, j) !== getEdgeIndex(j, i)`
+     */
+    getDirectedEdgeIndex(i: number, j: number): number,
     getVertexEdgeCount(i: number): number
 }
 
@@ -50,6 +59,13 @@ namespace IntAdjacencyGraph {
             return -1;
         }
 
+        getDirectedEdgeIndex(i: number, j: number): number {
+            for (let t = this.offset[i], _t = this.offset[i + 1]; t < _t; t++) {
+                if (this.b[t] === j) return t;
+            }
+            return -1;
+        }
+
         getVertexEdgeCount(i: number): number {
             return this.offset[i + 1] - this.offset[i];
         }

src/mol-view/stage.ts

@@ -11,15 +11,7 @@ import { MmcifUrlToModel, ModelToStructure, StructureToSpacefill, StructureToBal
 import { UrlEntity } from './state/entity';
 import { SpacefillProps } from 'mol-geo/representation/structure/spacefill';
 import { Context } from 'mol-app/context/context';
-
-// export const ColorTheme = {
-//     'atom-index': {},
-//     'chain-id': {},
-//     'element-symbol': {},
-//     'instance-index': {},
-//     'uniform': {}
-// }
-// export type ColorTheme = keyof typeof ColorTheme
+import { BallAndStickProps } from 'mol-geo/representation/structure/ball-and-stick';
 
 const spacefillProps: SpacefillProps = {
     doubleSided: true,
@@ -27,6 +19,14 @@ const spacefillProps: SpacefillProps = {
     colorTheme: { name: 'atom-index' }
 }
 
+const ballAndStickProps: BallAndStickProps = {
+    doubleSided: true,
+    detail: 1,
+    radialSegments: 8,
+    colorTheme: { name: 'chain-id' },
+    sizeTheme: { name: 'uniform', value: 0.25 },
+}
+
 export class Stage {
     viewer: Viewer
     ctx = new StateContext(Progress.format)
@@ -39,22 +39,27 @@ export class Stage {
         this.viewer = Viewer.create(canvas, container)
         this.viewer.animate()
         this.ctx.viewer = this.viewer
-        //this.loadPdbid('1jj2')
-        this.loadMmcifUrl(`../../examples/1cbs_full.bcif`)
+
+        // this.loadPdbid('1jj2')
+        this.loadPdbid('4umt') // ligand has bond with order 3
+        // this.loadPdbid('1crn')
+        // this.loadMmcifUrl(`../../examples/1cbs_full.bcif`)
     }
 
     async loadMmcifUrl (url: string) {
         const urlEntity = UrlEntity.ofUrl(this.ctx, url)
         const modelEntity = await MmcifUrlToModel.apply(this.ctx, urlEntity)
         const structureEntity = await ModelToStructure.apply(this.ctx, modelEntity)
+
         StructureToSpacefill.apply(this.ctx, structureEntity, { ...spacefillProps, visible: false })
-        StructureToBallAndStick.apply(this.ctx, structureEntity, spacefillProps) // TODO props
+        StructureToBallAndStick.apply(this.ctx, structureEntity, ballAndStickProps)
 
         this.globalContext.components.sequenceView.setState({ structure: structureEntity.value });
     }
 
     loadPdbid (pdbid: string) {
-        return this.loadMmcifUrl(`https://files.rcsb.org/download/${pdbid}.cif`)
+        return this.loadMmcifUrl(`http://www.ebi.ac.uk/pdbe/static/entry/${pdbid}_updated.cif`)
+        // return this.loadMmcifUrl(`https://files.rcsb.org/download/${pdbid}.cif`)
     }
 
     dispose () {