wip, cartoon

src/mol-geo/representation/structure/visual/cross-link-restraint-cylinder.ts

@@ -8,7 +8,7 @@ import { ValueCell } from 'mol-util/value-cell'
 
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object'
 import { Link, Structure } from 'mol-model/structure';
-import { DefaultStructureProps, StructureVisual } from '../index';
+import { DefaultStructureProps, StructureVisual } from '..';
 import { RuntimeContext } from 'mol-task'
 import { LinkCylinderProps, DefaultLinkCylinderProps, createLinkCylinderMesh } from './util/link';
 import { MeshValues } from 'mol-gl/renderable';

src/mol-geo/representation/structure/visual/element-point.ts

@@ -11,7 +11,7 @@ import { Unit } from 'mol-model/structure';
 import { RuntimeContext } from 'mol-task'
 import { fillSerial } from 'mol-gl/renderable/util';
 
-import { UnitsVisual, DefaultStructureProps } from '../index';
+import { UnitsVisual, DefaultStructureProps } from '..';
 import VertexMap from '../../../shape/vertex-map';
 import { SizeTheme } from '../../../theme';
 import { markElement, getElementLoci } from './util/element';

src/mol-geo/representation/structure/visual/element-sphere.ts

@@ -9,7 +9,7 @@ import { ValueCell } from 'mol-util/value-cell'
 
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object'
 import { Unit } from 'mol-model/structure';
-import { DefaultStructureProps, UnitsVisual } from '../index';
+import { DefaultStructureProps, UnitsVisual } from '..';
 import { RuntimeContext } from 'mol-task'
 import { createTransforms, createColors } from './util/common';
 import { createElementSphereMesh, markElement, getElementRadius, getElementLoci } from './util/element';

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

@@ -8,7 +8,7 @@ import { ValueCell } from 'mol-util/value-cell'
 
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object'
 import { Link, Structure } from 'mol-model/structure';
-import { DefaultStructureProps, StructureVisual } from '../index';
+import { DefaultStructureProps, StructureVisual } from '..';
 import { RuntimeContext } from 'mol-task'
 import { LinkCylinderProps, DefaultLinkCylinderProps, createLinkCylinderMesh } from './util/link';
 import { MeshValues } from 'mol-gl/renderable';

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

@@ -9,7 +9,7 @@ import { ValueCell } from 'mol-util/value-cell'
 
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object'
 import { Unit, Link } from 'mol-model/structure';
-import { UnitsVisual, DefaultStructureProps } from '../index';
+import { UnitsVisual, DefaultStructureProps } from '..';
 import { RuntimeContext } from 'mol-task'
 import { DefaultLinkCylinderProps, LinkCylinderProps, createLinkCylinderMesh } from './util/link';
 import { MeshValues } from 'mol-gl/renderable';

src/mol-geo/representation/structure/visual/polymer-backbone-cylinder.ts

@@ -8,7 +8,7 @@ import { ValueCell } from 'mol-util/value-cell'
 
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object'
 import { Unit } from 'mol-model/structure';
-import { DefaultStructureProps, UnitsVisual } from '../index';
+import { DefaultStructureProps, UnitsVisual } from '..';
 import { RuntimeContext } from 'mol-task'
 import { createTransforms, createColors } from './util/common';
 import { deepEqual } from 'mol-util';
@@ -26,8 +26,8 @@ import { getElementLoci, markElement } from './util/element';
 
 async function createPolymerBackboneCylinderMesh(ctx: RuntimeContext, unit: Unit, mesh?: Mesh) {
     const polymerElementCount = getPolymerElementCount(unit)
-    console.log('polymerElementCount backbone', polymerElementCount)
     if (!polymerElementCount) return Mesh.createEmpty(mesh)
+    console.log('polymerElementCount backbone', polymerElementCount)
 
     // TODO better vertex count estimates
     const builder = MeshBuilder.create(polymerElementCount * 30, polymerElementCount * 30 / 2, mesh)
@@ -80,7 +80,7 @@ export function PolymerBackboneVisual(): UnitsVisual<PolymerBackboneProps> {
             mesh = unitKinds.includes(unit.kind)
                 ? await createPolymerBackboneCylinderMesh(ctx, unit, mesh)
                 : Mesh.createEmpty(mesh)
-            console.log(mesh)
+            // console.log(mesh)
 
             const transforms = createTransforms(group)
             const color = createColors(group, elementCount, colorTheme)

src/mol-geo/representation/structure/visual/polymer-trace-mesh.ts

@@ -8,7 +8,7 @@ import { ValueCell } from 'mol-util/value-cell'
 
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object'
 import { Unit, Element } from 'mol-model/structure';
-import { DefaultStructureProps, UnitsVisual } from '../index';
+import { DefaultStructureProps, UnitsVisual } from '..';
 import { RuntimeContext } from 'mol-task'
 import { createTransforms, createColors } from './util/common';
 import { markElement } from './util/element';
@@ -25,11 +25,58 @@ import { createMeshValues, updateMeshValues, updateRenderableState, createRender
 import { MeshBuilder } from '../../../shape/mesh-builder';
 import { getPolymerElementCount, PolymerTraceIterator } from './util/polymer';
 import { Vec3 } from 'mol-math/linear-algebra';
+import { SecondaryStructureType } from 'mol-model/structure/model/types';
+// import { radToDeg } from 'mol-math/misc';
+
+const tmpNormal = Vec3.zero()
+const tangentVec = Vec3.zero()
+const normalVec = Vec3.zero()
+const binormalVec = Vec3.zero()
+
+const prevNormal = Vec3.zero()
+
+const orthogonalizeTmpVec = Vec3.zero()
+/** Get a vector that is similar to b but orthogonal to a */
+function orthogonalize(out: Vec3, a: Vec3, b: Vec3) {
+    Vec3.normalize(orthogonalizeTmpVec, Vec3.cross(orthogonalizeTmpVec, a, b))
+    Vec3.normalize(out, Vec3.cross(out, orthogonalizeTmpVec, a))
+    return out
+}
+
+function interpolateNormals(controlPoints: Helpers.NumberArray, tangentVectors: Helpers.NumberArray, normalVectors: Helpers.NumberArray, binormalVectors: Helpers.NumberArray, firstNormalVector: Vec3, lastNormalVector: Vec3) {
+    const n = controlPoints.length / 3
+    // const n1 = n - 1
+
+    // const angle = radToDeg(Math.acos(Vec3.dot(firstNormalVector, lastNormalVector)))
+    // console.log('angle', angle)
+    if (Vec3.dot(firstNormalVector, lastNormalVector) < 0) {
+        Vec3.scale(lastNormalVector, lastNormalVector, -1)
+        // console.log('flipped last normal vector')
+    }
+
+    Vec3.copy(prevNormal, firstNormalVector)
+
+    for (let i = 0; i < n; ++i) {
+        const t = i === 0 ? 0 : 1 / (n - i)
+        Vec3.normalize(tmpNormal, Vec3.slerp(tmpNormal, prevNormal, lastNormalVector, t))
+
+        Vec3.fromArray(tangentVec, tangentVectors, i * 3)
 
-export function reflect(target: Vec3, p1: Vec3, p2: Vec3, amount: number) {
-    target[0] = p1[0] - amount * (p2[0] - p1[0])
-    target[1] = p1[1] - amount * (p2[1] - p1[1])
-    target[2] = p1[2] - amount * (p2[2] - p1[2])
+        orthogonalize(normalVec, tangentVec, tmpNormal)
+        Vec3.toArray(normalVec, normalVectors, i * 3)
+
+        // const deltaAngle = radToDeg(Math.acos(Vec3.dot(prevNormal, normalVec)))
+        // if (deltaAngle > (angle / n1) * 5 && deltaAngle > 20) {
+        //     console.warn(i, 'large delta angle', deltaAngle)
+        // }
+        // if (Vec3.dot(normalVec, prevNormal) < 0) {
+        //     console.warn(i, 'flip compared to prev', radToDeg(Math.acos(Vec3.dot(prevNormal, normalVec))))
+        // }
+        Vec3.copy(prevNormal, normalVec)
+        
+        Vec3.normalize(binormalVec, Vec3.cross(binormalVec, tangentVec, normalVec))
+        Vec3.toArray(binormalVec, binormalVectors, i * 3)
+    }
 }
 
 async function createPolymerTraceMesh(ctx: RuntimeContext, unit: Unit, mesh?: Mesh) {
@@ -39,26 +86,23 @@ async function createPolymerTraceMesh(ctx: RuntimeContext, unit: Unit, mesh?: Me
 
     // TODO better vertex count estimates
     const builder = MeshBuilder.create(polymerElementCount * 30, polymerElementCount * 30 / 2, mesh)
-    const linearSegments = 5
-    const radialSegments = 8
+    const linearSegments = 12
+    const radialSegments = 16
     const tension = 0.9
 
-    const tA = Vec3.zero()
+    const tanA = Vec3.zero()
+    const tanB = Vec3.zero()
+
     const tB = Vec3.zero()
-    const dA = Vec3.zero()
-    const dB = Vec3.zero()
-    const torsionVec = Vec3.zero()
-    const initialTorsionVec = Vec3.zero()
     const tangentVec = Vec3.zero()
-    const normalVec = Vec3.zero()
 
     const tmp = Vec3.zero()
-    const reflectedControlPoint = Vec3.zero()
 
     const pn = (linearSegments + 1) * 3
     const controlPoints = new Float32Array(pn)
-    const torsionVectors = new Float32Array(pn)
+    const tangentVectors = new Float32Array(pn)
     const normalVectors = new Float32Array(pn)
+    const binormalVectors = new Float32Array(pn)
 
     let i = 0
     const polymerTraceIt = PolymerTraceIterator(unit)
@@ -66,120 +110,88 @@ async function createPolymerTraceMesh(ctx: RuntimeContext, unit: Unit, mesh?: Me
         const v = polymerTraceIt.move()
         builder.setId(v.index)
 
-        Vec3.spline(tB, v.t1, v.t2, v.t3, v.t4, 0.5, tension)
-        Vec3.spline(dA, v.d12, v.d23, v.d34, v.d45, 0.5, tension)
-
-        Vec3.normalize(initialTorsionVec, Vec3.sub(initialTorsionVec, tB, dB))
-
-        Vec3.toArray(tB, controlPoints, 0)
-        Vec3.normalize(torsionVec, Vec3.sub(torsionVec, tB, dB))
-        Vec3.toArray(torsionVec, torsionVectors, 0)
-        // approximate tangent as direction to previous control point
-        Vec3.normalize(tangentVec, Vec3.sub(tangentVec, tB, tA))
-        Vec3.normalize(normalVec, Vec3.cross(normalVec, tangentVec, torsionVec))
-        Vec3.toArray(normalVec, normalVectors, 0)
-
-        //
-
-        const t12 = Vec3.zero()
-        const t23 = Vec3.zero()
-        const t34 = Vec3.zero()
-        const t45 = Vec3.zero()
-        Vec3.spline(t12, v.t0, v.t1, v.t2, v.t3, 0.5, tension)
-        Vec3.spline(t23, v.t1, v.t2, v.t3, v.t4, 0.5, tension)
-        Vec3.spline(t34, v.t2, v.t3, v.t4, v.t5, 0.5, tension)
-        Vec3.spline(t45, v.t3, v.t4, v.t5, v.t6, 0.5, tension)
-
-        // const dp12 = Vec3.zero()
-        // const dp23 = Vec3.zero()
-        // const dp34 = Vec3.zero()
-        // const dp45 = Vec3.zero()
-        // Vec3.projectPointOnVector(dp12, v.d12, t12, v.t1)
-        // Vec3.projectPointOnVector(dp23, v.d23, t23, v.t2)
-        // Vec3.projectPointOnVector(dp34, v.d34, t34, v.t3)
-        // Vec3.projectPointOnVector(dp45, v.d45, t45, v.t4)
-
-        const td12 = Vec3.zero()
-        const td23 = Vec3.zero()
-        const td34 = Vec3.zero()
-        const td45 = Vec3.zero()
-        Vec3.normalize(td12, Vec3.sub(td12, t12, v.d12))
-        Vec3.scaleAndAdd(v.d12, t12, td12, 1)
-        Vec3.normalize(td23, Vec3.sub(td23, t23, v.d23))
-        if (Vec3.dot(td12, td23) < 0) {
-            Vec3.scaleAndAdd(v.d23, t23, td23, -1)
-            console.log('foo td0 td1')
-        } else {
-            Vec3.scaleAndAdd(v.d23, t23, td23, 1)
-        }
-        Vec3.normalize(td34, Vec3.sub(td34, t34, v.d34))
-        if (Vec3.dot(td12, td34) < 0) {
-            Vec3.scaleAndAdd(v.d34, t34, td34, -1)
-            console.log('foo td1 td2')
-        } else {
-            Vec3.scaleAndAdd(v.d34, t34, td34, 1)
-        }
-        Vec3.normalize(td45, Vec3.sub(td45, t45, v.d45))
-        if (Vec3.dot(td12, td45) < 0) {
-            Vec3.scaleAndAdd(v.d45, t45, td45, -1)
-            console.log('foo td2 td3')
-        } else {
-            Vec3.scaleAndAdd(v.d45, t45, td45, 1)
-        }
-
-        // console.log(td0, td1, td2, td3)
-
-        builder.addIcosahedron(t12, 0.3, 1)
-        builder.addIcosahedron(t23, 0.3, 1)
-        builder.addIcosahedron(t34, 0.3, 1)
-        builder.addIcosahedron(t45, 0.3, 1)
-
-        // builder.addIcosahedron(dp12, 0.3, 1)
-        // builder.addIcosahedron(dp23, 0.3, 1)
-        // builder.addIcosahedron(dp34, 0.3, 1)
-        // builder.addIcosahedron(dp45, 0.3, 1)
-
-        builder.addIcosahedron(v.d12, 0.3, 1)
-        builder.addIcosahedron(v.d23, 0.3, 1)
-        builder.addIcosahedron(v.d34, 0.3, 1)
-        builder.addIcosahedron(v.d45, 0.3, 1)
-
-        for (let j = 1; j <= linearSegments; ++j) {
+        for (let j = 0; j <= linearSegments; ++j) {
             const t = j * 1.0 / linearSegments;
-            Vec3.copy(tA, tB)
             // if ((v.last && t > 0.5) || (v.first && t < 0.5)) break
 
             if (t < 0.5) {
-                Vec3.spline(tB, v.t1, v.t2, v.t3, v.t4, t + 0.5, tension)
+                Vec3.spline(tB, v.t0, v.t1, v.t2, v.t3, t + 0.5, tension)
+                Vec3.spline(tanA, v.t0, v.t1, v.t2, v.t3, t + 0.5 + 0.01, tension)
+                Vec3.spline(tanB, v.t0, v.t1, v.t2, v.t3, t + 0.5 - 0.01, tension)
             } else {
-                Vec3.spline(tB, v.t2, v.t3, v.t4, v.t5, t - 0.5, tension)
+                Vec3.spline(tB, v.t1, v.t2, v.t3, v.t4, t - 0.5, tension)
+                Vec3.spline(tanA, v.t1, v.t2, v.t3, v.t4, t - 0.5 + 0.01, tension)
+                Vec3.spline(tanB, v.t1, v.t2, v.t3, v.t4, t - 0.5 - 0.01, tension)
             }
-            Vec3.spline(dB, v.d12, v.d23, v.d34, v.d45, t, tension)
-
-            // reflect(reflectedControlPoint, tB, tA, 1)
             Vec3.toArray(tB, controlPoints, j * 3)
+            Vec3.normalize(tangentVec, Vec3.sub(tangentVec, tanA, tanB))
+            Vec3.toArray(tangentVec, tangentVectors, j * 3)
+        }
 
-            Vec3.normalize(torsionVec, Vec3.sub(torsionVec, tB, dB))
-            // if (Vec3.dot(initialTorsionVec, torsionVec) < 0) Vec3.scale(torsionVec, torsionVec, -1)
-            Vec3.toArray(torsionVec, torsionVectors, j * 3)
-
-            // approximate tangent as direction to previous control point
-            Vec3.normalize(tangentVec, Vec3.sub(tangentVec, tB, tA))
-            Vec3.normalize(normalVec, Vec3.cross(normalVec, tangentVec, torsionVec))
-            Vec3.toArray(normalVec, normalVectors, j * 3)
-
-            // TODO size theme
-            // builder.addCylinder(tA, tB, 1.0, { radiusTop: 0.3, radiusBottom: 0.3 })
-
-            builder.addIcosahedron(dB, 0.1, 1)
-
-            builder.addCylinder(tB, Vec3.add(tmp, tB, torsionVec), 1.0, { radiusTop: 0.1, radiusBottom: 0.1 })
-            // builder.addCylinder(tB, Vec3.add(tmp, tB, normalVec), 1.0, { radiusTop: 0.1, radiusBottom: 0.1 })
-
-            console.log(tA, tB)
+        const firstControlPoint = Vec3.zero()
+        const lastControlPoint = Vec3.zero()
+        const firstTangentVec = Vec3.zero()
+        const lastTangentVec = Vec3.zero()
+        const firstNormalVec = Vec3.zero()
+        const lastNormalVec = Vec3.zero()
+        const firstDirPoint = Vec3.zero()
+        const lastDirPoint = Vec3.zero()
+
+        Vec3.fromArray(firstControlPoint, controlPoints, 0)
+        Vec3.fromArray(lastControlPoint, controlPoints, linearSegments * 3)
+        Vec3.fromArray(firstTangentVec, tangentVectors, 0)
+        Vec3.fromArray(lastTangentVec, tangentVectors, linearSegments * 3)
+        Vec3.copy(firstDirPoint, v.d12)
+        Vec3.copy(lastDirPoint, v.d23)
+
+        Vec3.normalize(tmpNormal, Vec3.sub(tmp, firstControlPoint, firstDirPoint))
+        orthogonalize(firstNormalVec, firstTangentVec, tmpNormal)
+
+        Vec3.normalize(tmpNormal, Vec3.sub(tmp, lastControlPoint, lastDirPoint))
+        orthogonalize(lastNormalVec, lastTangentVec, tmpNormal)
+
+        // console.log('ELEMENT', i)
+        interpolateNormals(controlPoints, tangentVectors, normalVectors, binormalVectors, firstNormalVec, lastNormalVec)
+
+        // const controlPoint = Vec3.zero()
+        // for (let j = 0; j <= linearSegments; ++j) {
+        //     Vec3.fromArray(controlPoint, controlPoints, j * 3)
+        //     Vec3.fromArray(normalVec, normalVectors, j * 3)
+        //     Vec3.fromArray(binormalVec, binormalVectors, j * 3)
+        //     Vec3.fromArray(tangentVec, tangentVectors, j * 3)
+        //     builder.addIcosahedron(controlPoint, 0.1, 1)
+        //     builder.addCylinder(
+        //         controlPoint,
+        //         Vec3.add(tmp, controlPoint, normalVec),
+        //         1.5,
+        //         { radiusTop: 0.07, radiusBottom: 0.07 }
+        //     )
+        //     builder.addCylinder(
+        //         controlPoint,
+        //         Vec3.add(tmp, controlPoint, binormalVec),
+        //         0.8,
+        //         { radiusTop: 0.07, radiusBottom: 0.07 }
+        //     )
+        //     builder.addCylinder(
+        //         controlPoint,
+        //         Vec3.add(tmp, controlPoint, tangentVec),
+        //         j === 0 ? 2 : 1.5,
+        //         { radiusTop: 0.03, radiusBottom: 0.03 }
+        //     )
+        // }
+
+        let width = 0.2
+        let height = 0.2
+        if (
+            SecondaryStructureType.is(v.secStrucType, SecondaryStructureType.Flag.Beta) ||
+            SecondaryStructureType.is(v.secStrucType, SecondaryStructureType.Flag.Helix)
+        ) {
+            width = 0.2
+            height = 0.6
         }
 
-        builder.addTube(controlPoints, torsionVectors, normalVectors, linearSegments, radialSegments)
+        // TODO size theme
+        builder.addTube(controlPoints, normalVectors, binormalVectors, linearSegments, radialSegments, width, height)
 
         if (i % 10000 === 0 && ctx.shouldUpdate) {
             await ctx.update({ message: 'Polymer trace mesh', current: i, max: polymerElementCount });

src/mol-geo/representation/structure/visual/util/polymer.ts

@@ -5,45 +5,46 @@
  */
 
 import { Unit, Element, StructureProperties, Model } from 'mol-model/structure';
-import { Segmentation, Interval } from 'mol-data/int';
-import { MoleculeType } from 'mol-model/structure/model/types';
+import { Segmentation, OrderedSet, Interval } from 'mol-data/int';
+import { MoleculeType, SecondaryStructureType } from 'mol-model/structure/model/types';
 import Iterator from 'mol-data/iterator';
-import { SegmentIterator } from 'mol-data/int/impl/segmentation';
 import { Vec3 } from 'mol-math/linear-algebra';
 import { SymmetryOperator } from 'mol-math/geometry';
+import SortedRanges from 'mol-data/int/sorted-ranges';
 
-// type TraceMap = Map<number, number>
-
-// interface TraceMaps {
-//     atomic: TraceMap
-//     spheres: TraceMap
-//     gaussians: TraceMap
-// }
-
-// function calculateTraceMaps (model: Model): TraceMaps {
-
-// }
+export function getPolymerRanges(unit: Unit) {
+    switch (unit.kind) {
+        case Unit.Kind.Atomic: return unit.model.atomicHierarchy.polymerRanges
+        case Unit.Kind.Spheres: return unit.model.coarseHierarchy.spheres.polymerRanges
+        case Unit.Kind.Gaussians: return unit.model.coarseHierarchy.gaussians.polymerRanges
+    }
+}
 
 export function getPolymerElementCount(unit: Unit) {
     let count = 0
     const { elements } = unit
-    const l = Element.Location(unit)
-    if (Unit.isAtomic(unit)) {
-        const { polymerSegments, residueSegments } = unit.model.atomicHierarchy
-        const polymerIt = Segmentation.transientSegments(polymerSegments, elements);
-        const residuesIt = Segmentation.transientSegments(residueSegments, elements);
-        while (polymerIt.hasNext) {
-            residuesIt.setSegment(polymerIt.move());
-            while (residuesIt.hasNext) {
-                residuesIt.move();
-                count++
+    const polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), elements)
+    switch (unit.kind) {
+        case Unit.Kind.Atomic:
+            const residueIt = Segmentation.transientSegments(unit.model.atomicHierarchy.residueSegments, elements)
+            while (polymerIt.hasNext) {
+                const polymerSegment = polymerIt.move()
+                residueIt.setSegment(polymerSegment)
+                while (residueIt.hasNext) {
+                    const residueSegment = residueIt.move()
+                    const { start, end } = residueSegment
+                    if (OrderedSet.areIntersecting(Interval.ofBounds(elements[start], elements[end - 1]), elements)) ++count
+                }
             }
-        }
-    } else if (Unit.isSpheres(unit)) {
-        for (let i = 0, il = elements.length; i < il; ++i) {
-            l.element = elements[i]
-            if (StructureProperties.entity.type(l) === 'polymer') count++
-        }
+            break
+        case Unit.Kind.Spheres:
+        case Unit.Kind.Gaussians:
+            while (polymerIt.hasNext) {
+                const { start, end } = polymerIt.move()
+                // TODO add OrderedSet.intersectionSize
+                count += OrderedSet.size(OrderedSet.intersect(Interval.ofBounds(elements[start], elements[end - 1]), elements))
+            }
+            break
     }
     return count
 }
@@ -99,7 +100,7 @@ function getTraceElement2(model: Model, residueModelSegment: Segmentation.Segmen
     for (let j = residueModelSegment.start, _j = residueModelSegment.end; j < _j; j++) {
         if (model.atomicHierarchy.atoms.label_atom_id.value(j) === traceName) return j
     }
-    console.log('trace name element not found', { ...residueModelSegment })
+    console.log(`trace name element "${traceName}" not found`, { ...residueModelSegment })
     return residueModelSegment.start
 }
 
@@ -115,22 +116,31 @@ function getDirectionName2(model: Model, residueModelIndex: number) {
         traceName = 'O4\''
     } else if (moleculeType === MoleculeType.RNA) {
         traceName = 'C3\''
+    } else {
+        console.log('molecule type unknown', Number.isInteger(residueModelIndex) ? compId : residueModelIndex, chemCompMap)
     }
     return traceName
 }
 
+// function residueLabel(model: Model, rI: number) {
+//     const { residues, chains, residueSegments, chainSegments } = model.atomicHierarchy
+//     const { label_comp_id, label_seq_id } = residues
+//     const { label_asym_id } = chains
+//     const cI = chainSegments.segmentMap[residueSegments.segments[rI]]
+//     return `${label_asym_id.value(cI)} ${label_comp_id.value(rI)} ${label_seq_id.value(rI)}`
+// }
+
 function getDirectionElement2(model: Model, residueModelSegment: Segmentation.Segment<Element>) {
     const traceName = getDirectionName2(model, residueModelSegment.index)
 
     for (let j = residueModelSegment.start, _j = residueModelSegment.end; j < _j; j++) {
         if (model.atomicHierarchy.atoms.label_atom_id.value(j) === traceName) return j
     }
-    console.log('direction name element not found', { ...residueModelSegment })
+    // console.log('direction name element not found', { ...residueModelSegment })
     return residueModelSegment.start
 }
 
 
-
 /** Iterates over consecutive pairs of residues/coarse elements in polymers */
 export function PolymerBackboneIterator(unit: Unit): Iterator<PolymerBackbonePair> {
     switch (unit.kind) {
@@ -166,8 +176,8 @@ const enum AtomicPolymerBackboneIteratorState { nextPolymer, firstResidue, nextR
 export class AtomicPolymerBackboneIterator<T extends number = number> implements Iterator<PolymerBackbonePair> {
     private value: PolymerBackbonePair
 
-    private polymerIt: SegmentIterator<Element>
-    private residueIt: SegmentIterator<Element>
+    private polymerIt: SortedRanges.Iterator<Element>
+    private residueIt: Segmentation.Iterator<Element>
     private polymerSegment: Segmentation.Segment<Element>
     private state: AtomicPolymerBackboneIteratorState = AtomicPolymerBackboneIteratorState.nextPolymer
     private pos: SymmetryOperator.CoordinateMapper
@@ -178,26 +188,25 @@ export class AtomicPolymerBackboneIterator<T extends number = number> implements
         const { residueIt, polymerIt, value, pos } = this
 
         if (this.state === AtomicPolymerBackboneIteratorState.nextPolymer) {
-            if (polymerIt.hasNext) {
+            while (polymerIt.hasNext) {
                 this.polymerSegment = polymerIt.move();
+                // console.log('polymerSegment', this.polymerSegment)
                 residueIt.setSegment(this.polymerSegment);
-                this.state = AtomicPolymerBackboneIteratorState.firstResidue
-            }
-        }
 
-        if (this.state === AtomicPolymerBackboneIteratorState.firstResidue) {
-            const residueSegment = residueIt.move();
-            if (residueIt.hasNext) {
-                value.indexB = setTraceElement(value.centerB, residueSegment)
-                pos(value.centerB.element, value.posB)
-                this.state = AtomicPolymerBackboneIteratorState.nextResidue
-            } else {
-                this.state = AtomicPolymerBackboneIteratorState.nextPolymer
+                const residueSegment = residueIt.move();
+                // console.log('first residueSegment', residueSegment, residueIt.hasNext)
+                if (residueIt.hasNext) {
+                    value.indexB = setTraceElement(value.centerB, residueSegment)
+                    pos(value.centerB.element, value.posB)
+                    this.state = AtomicPolymerBackboneIteratorState.nextResidue
+                    break
+                }
             }
         }
 
         if (this.state === AtomicPolymerBackboneIteratorState.nextResidue) {
             const residueSegment = residueIt.move();
+            // console.log('next residueSegment', residueSegment)
             value.centerA.element = value.centerB.element
             value.indexA = value.indexB
             Vec3.copy(value.posA, value.posB)
@@ -205,31 +214,37 @@ export class AtomicPolymerBackboneIterator<T extends number = number> implements
             pos(value.centerB.element, value.posB)
 
             if (!residueIt.hasNext) {
+                // TODO need to advance to a polymer that has two or more residues (can't assume it has)
                 this.state = AtomicPolymerBackboneIteratorState.nextPolymer
             }
         }
 
         this.hasNext = residueIt.hasNext || polymerIt.hasNext
+        
+        // console.log('hasNext', this.hasNext)
+        // console.log('value', this.value)
 
         return this.value;
     }
 
     constructor(unit: Unit.Atomic) {
-        const { polymerSegments, residueSegments } = unit.model.atomicHierarchy
-        this.polymerIt = Segmentation.transientSegments(polymerSegments, unit.elements);
-        this.residueIt = Segmentation.transientSegments(residueSegments, unit.elements);
+        const { residueSegments } = unit.model.atomicHierarchy
+        // console.log('unit.elements', OrderedSet.toArray(unit.elements))
+        this.polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), unit.elements)
+        this.residueIt = Segmentation.transientSegments(residueSegments, unit.elements)
         this.pos = unit.conformation.invariantPosition
         this.value = createPolymerBackbonePair(unit)
-        this.hasNext = this.residueIt.hasNext || this.polymerIt.hasNext
+
+        this.hasNext = this.residueIt.hasNext && this.polymerIt.hasNext
     }
 }
 
-const enum CoarsePolymerBackboneIteratorState { nextPolymer, firstElement, nextElement }
+const enum CoarsePolymerBackboneIteratorState { nextPolymer, nextElement }
 
 export class CoarsePolymerBackboneIterator<T extends number = number> implements Iterator<PolymerBackbonePair> {
     private value: PolymerBackbonePair
 
-    private polymerIt: SegmentIterator<Element>
+    private polymerIt: SortedRanges.Iterator<Element>
     private polymerSegment: Segmentation.Segment<Element>
     private state: CoarsePolymerBackboneIteratorState = CoarsePolymerBackboneIteratorState.nextPolymer
     private pos: SymmetryOperator.CoordinateMapper
@@ -243,25 +258,21 @@ export class CoarsePolymerBackboneIterator<T extends number = number> implements
         if (this.state === CoarsePolymerBackboneIteratorState.nextPolymer) {
             if (polymerIt.hasNext) {
                 this.polymerSegment = polymerIt.move();
+                console.log('polymer', this.polymerSegment)
                 this.elementIndex = this.polymerSegment.start
-                this.state = CoarsePolymerBackboneIteratorState.firstElement
+                this.elementIndex += 1
+                if (this.elementIndex + 1 < this.polymerSegment.end) {
+                    value.centerB.element = value.centerB.unit.elements[this.elementIndex]
+                    value.indexB = this.elementIndex
+                    pos(value.centerB.element, value.posB)
+
+                    this.state = CoarsePolymerBackboneIteratorState.nextElement
+                } else {
+                    this.state = CoarsePolymerBackboneIteratorState.nextPolymer
+                }
             }
         }
 
-        if (this.state === CoarsePolymerBackboneIteratorState.firstElement) {
-            this.elementIndex += 1
-            if (this.elementIndex + 1 < this.polymerSegment.end) {
-                value.centerB.element = value.centerB.unit.elements[this.elementIndex]
-                value.indexB = this.elementIndex
-                pos(value.centerB.element, value.posB)
-
-                this.state = CoarsePolymerBackboneIteratorState.nextElement
-            } else {
-                this.state = CoarsePolymerBackboneIteratorState.nextPolymer
-            }
-
-        }
-
         if (this.state === CoarsePolymerBackboneIteratorState.nextElement) {
             this.elementIndex += 1
             value.centerA.element = value.centerB.element
@@ -282,14 +293,13 @@ export class CoarsePolymerBackboneIterator<T extends number = number> implements
     }
 
     constructor(unit: Unit.Spheres | Unit.Gaussians) {
-        const { polymerSegments } = Unit.isSpheres(unit)
-            ? unit.model.coarseHierarchy.spheres
-            : unit.model.coarseHierarchy.gaussians
-        this.polymerIt = Segmentation.transientSegments(polymerSegments, unit.elements);
+        this.polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), unit.elements);
 
         this.pos = unit.conformation.invariantPosition
         this.value = createPolymerBackbonePair(unit)
 
+        console.log('CoarsePolymerBackboneIterator', this.polymerIt.hasNext)
+
         this.hasNext = this.polymerIt.hasNext
     }
 }
@@ -315,19 +325,16 @@ interface PolymerTraceElement {
     index: number
     first: boolean
     last: boolean
+    secStrucType: SecondaryStructureType
 
     t0: Vec3
     t1: Vec3
     t2: Vec3
     t3: Vec3
     t4: Vec3
-    t5: Vec3
-    t6: Vec3
 
     d12: Vec3
     d23: Vec3
-    d34: Vec3
-    d45: Vec3
 }
 
 function createPolymerTraceElement (unit: Unit): PolymerTraceElement {
@@ -336,19 +343,16 @@ function createPolymerTraceElement (unit: Unit): PolymerTraceElement {
         index: 0,
         first: false,
         last: false,
+        secStrucType: SecondaryStructureType.create(SecondaryStructureType.Flag.NA),
 
         t0: Vec3.zero(),
         t1: Vec3.zero(),
         t2: Vec3.zero(),
         t3: Vec3.zero(),
         t4: Vec3.zero(),
-        t5: Vec3.zero(),
-        t6: Vec3.zero(),
 
         d12: Vec3.zero(),
         d23: Vec3.zero(),
-        d34: Vec3.zero(),
-        d45: Vec3.zero(),
     }
 }
 
@@ -368,8 +372,8 @@ function setSegment (outSegment: Segmentation.Segment<Element>, index: number, s
 export class AtomicPolymerTraceIterator<T extends number = number> implements Iterator<PolymerTraceElement> {
     private value: PolymerTraceElement
 
-    private polymerIt: SegmentIterator<Element>
-    private residueIt: SegmentIterator<Element>
+    private polymerIt: SortedRanges.Iterator<Element>
+    private residueIt: Segmentation.Iterator<Element>
     private residueSegmentMin: number
     private residueSegmentMax: number
     private state: AtomicPolymerTraceIteratorState = AtomicPolymerTraceIteratorState.nextPolymer
@@ -388,9 +392,9 @@ export class AtomicPolymerTraceIterator<T extends number = number> implements It
     }
 
     updateResidueSegmentRange(polymerSegment: Segmentation.Segment<Element>) {
-        const { polymerSegments, residueSegments } = this.unit.model.atomicHierarchy
-        const sMin = polymerSegments.segments[polymerSegment.index]
-        const sMax = polymerSegments.segments[polymerSegment.index + 1] - 1
+        const { polymerRanges, residueSegments } = this.unit.model.atomicHierarchy
+        const sMin = polymerRanges[polymerSegment.index * 2]
+        const sMax = polymerRanges[polymerSegment.index * 2 + 1]
         this.residueSegmentMin = residueSegments.segmentMap[sMin]
         this.residueSegmentMax = residueSegments.segmentMap[sMax]
     }
@@ -399,11 +403,15 @@ export class AtomicPolymerTraceIterator<T extends number = number> implements It
         const { residueIt, polymerIt, value } = this
 
         if (this.state === AtomicPolymerTraceIteratorState.nextPolymer) {
-            if (polymerIt.hasNext) {
+            while (polymerIt.hasNext) {
                 const polymerSegment = polymerIt.move();
+                // console.log('polymerSegment', {...polymerSegment})
                 residueIt.setSegment(polymerSegment);
                 this.updateResidueSegmentRange(polymerSegment)
-                this.state = AtomicPolymerTraceIteratorState.nextResidue
+                if (residueIt.hasNext) {
+                    this.state = AtomicPolymerTraceIteratorState.nextResidue
+                    break
+                }
             }
         }
 
@@ -411,32 +419,26 @@ export class AtomicPolymerTraceIterator<T extends number = number> implements It
             const { tmpSegment, residueSegments, residueSegmentMin, residueSegmentMax } = this
             const residueSegment = residueIt.move();
             const resSegIdx = residueSegment.index
+            // console.log(residueLabel(this.unit.model, resSegIdx), resSegIdx, this.unit.model.properties.secondaryStructure.type[resSegIdx])
             value.index = setTraceElement(value.center, residueSegment)
 
-            setSegment(tmpSegment, resSegIdx - 3, residueSegments, residueSegmentMin, residueSegmentMax)
+            setSegment(tmpSegment, resSegIdx - 2, residueSegments, residueSegmentMin, residueSegmentMax)
             this.pos(value.t0, getTraceElement2(this.unit.model, tmpSegment))
 
-            setSegment(tmpSegment, resSegIdx - 2, residueSegments, residueSegmentMin, residueSegmentMax)
+            setSegment(tmpSegment, resSegIdx - 1, residueSegments, residueSegmentMin, residueSegmentMax)
             this.pos(value.t1, getTraceElement2(this.unit.model, tmpSegment))
             this.pos(value.d12, getDirectionElement2(this.unit.model, tmpSegment))
 
-            setSegment(tmpSegment, resSegIdx - 1, residueSegments, residueSegmentMin, residueSegmentMax)
+            setSegment(tmpSegment, resSegIdx, residueSegments, residueSegmentMin, residueSegmentMax)
+            value.secStrucType = this.unit.model.properties.secondaryStructure.type[resSegIdx]
             this.pos(value.t2, getTraceElement2(this.unit.model, tmpSegment))
             this.pos(value.d23, getDirectionElement2(this.unit.model, tmpSegment))
 
-            setSegment(tmpSegment, resSegIdx, residueSegments, residueSegmentMin, residueSegmentMax)
-            this.pos(value.t3, getTraceElement2(this.unit.model, tmpSegment))
-            this.pos(value.d34, getDirectionElement2(this.unit.model, tmpSegment))
-
             setSegment(tmpSegment, resSegIdx + 1, residueSegments, residueSegmentMin, residueSegmentMax)
-            this.pos(value.t4, getTraceElement2(this.unit.model, tmpSegment))
-            this.pos(value.d45, getDirectionElement2(this.unit.model, tmpSegment))
+            this.pos(value.t3, getTraceElement2(this.unit.model, tmpSegment))
 
             setSegment(tmpSegment, resSegIdx + 2, residueSegments, residueSegmentMin, residueSegmentMax)
-            this.pos(value.t5, getTraceElement2(this.unit.model, tmpSegment))
-
-            setSegment(tmpSegment, resSegIdx + 3, residueSegments, residueSegmentMin, residueSegmentMax)
-            this.pos(value.t6, getTraceElement2(this.unit.model, tmpSegment))
+            this.pos(value.t4, getTraceElement2(this.unit.model, tmpSegment))
 
             value.first = resSegIdx === residueSegmentMin
             value.last = resSegIdx === residueSegmentMax
@@ -452,18 +454,16 @@ export class AtomicPolymerTraceIterator<T extends number = number> implements It
     }
 
     constructor(unit: Unit.Atomic) {
-        const { polymerSegments, residueSegments } = unit.model.atomicHierarchy
-        this.polymerIt = Segmentation.transientSegments(polymerSegments, unit.elements);
+        const { residueSegments } = unit.model.atomicHierarchy
+        this.polymerIt = SortedRanges.transientSegments(getPolymerRanges(unit), unit.elements)
         this.residueIt = Segmentation.transientSegments(residueSegments, unit.elements);
         this.residueSegments = residueSegments
         this.value = createPolymerTraceElement(unit)
-        this.hasNext = this.residueIt.hasNext || this.polymerIt.hasNext
+        this.hasNext = this.residueIt.hasNext && this.polymerIt.hasNext
 
         this.tmpSegment = { index: 0, start: 0 as Element, end: 0 as Element }
 
         this.unit = unit
-        console.log('model', unit.model)
-        console.log('unit', unit)
     }
 }
 

src/mol-geo/shape/mesh-builder.ts

@@ -27,7 +27,7 @@ export interface MeshBuilder {
     addDoubleCylinder(start: Vec3, end: Vec3, lengthScale: number, shift: Vec3, props: CylinderProps): void
     addFixedCountDashedCylinder(start: Vec3, end: Vec3, lengthScale: number, segmentCount: number, props: CylinderProps): void
     addIcosahedron(center: Vec3, radius: number, detail: number): void
-    addTube(controlPoints: Helpers.NumberArray, torsionVectors: Helpers.NumberArray, normalVectors: Helpers.NumberArray, linearSegments: number, radialSegments: number): void
+    addTube(controlPoints: Helpers.NumberArray, torsionVectors: Helpers.NumberArray, normalVectors: Helpers.NumberArray, linearSegments: number, radialSegments: number, width: number, height: number): void
     setId(id: number): void
     getMesh(): Mesh
 }
@@ -176,13 +176,13 @@ export namespace MeshBuilder {
                 setIcosahedronMat(tmpIcosahedronMat, center)
                 add(tmpIcosahedronMat, vertices, normals, indices)
             },
-            addTube: (controlPoints: Helpers.NumberArray, torsionVectors: Helpers.NumberArray, normalVectors: Helpers.NumberArray, linearSegments: number, radialSegments: number) => {
-                console.log(controlPoints, torsionVectors, normalVectors, linearSegments, radialSegments)
+            addTube: (controlPoints: Helpers.NumberArray, normalVectors: Helpers.NumberArray, binormalVectors: Helpers.NumberArray, linearSegments: number, radialSegments: number, width: number, height: number) => {
+                // console.log(controlPoints, normalVectors, binormalVectors, linearSegments, radialSegments)
 
-                const ico = getIcosahedron({ radius: 0.1, detail: 1 })
+                // const ico = getIcosahedron({ radius: 0.1, detail: 1 })
 
-                const radialVector = Vec3.zero()
                 const normalVector = Vec3.zero()
+                const binormalVector = Vec3.zero()
                 const tempPos = Vec3.zero()
                 const a = Vec3.zero()
                 const b = Vec3.zero()
@@ -190,16 +190,14 @@ export namespace MeshBuilder {
                 const v = Vec3.zero()
 
                 const waveFactor = 1
-                const width = 0.6
-                const height = 0.2
 
                 const vertexCount = vertices.elementCount
                 const di = 1 / linearSegments
 
                 for (let i = 0; i <= linearSegments; ++i) {
                     const i3 = i * 3
-                    Vec3.fromArray(u, torsionVectors, i3)
-                    Vec3.fromArray(v, normalVectors, i3)
+                    Vec3.fromArray(u, normalVectors, i3)
+                    Vec3.fromArray(v, binormalVectors, i3)
 
                     const tt = di * i - 0.5;
                     const ff = 1 + (waveFactor - 1) * (Math.cos(2 * Math.PI * tt) + 1);
@@ -211,7 +209,7 @@ export namespace MeshBuilder {
                         Vec3.copy(a, u)
                         Vec3.copy(b, v)
                         Vec3.add(
-                            radialVector,
+                            normalVector,
                             Vec3.scale(a, a, w * Math.cos(t)),
                             Vec3.scale(b, b, h * Math.sin(t))
                         )
@@ -219,14 +217,14 @@ export namespace MeshBuilder {
                         Vec3.copy(a, u)
                         Vec3.copy(b, v)
                         Vec3.add(
-                            normalVector,
+                            binormalVector,
                             Vec3.scale(a, a, h * Math.cos(t)),
                             Vec3.scale(b, b, w * Math.sin(t))
                         )
                         Vec3.normalize(normalVector, normalVector)
 
                         Vec3.fromArray(tempPos, controlPoints, i3)
-                        Vec3.add(tempPos, tempPos, radialVector)
+                        Vec3.add(tempPos, tempPos, binormalVector)
 
                         ChunkedArray.add3(vertices, tempPos[0], tempPos[1], tempPos[2]);
                         ChunkedArray.add3(normals, normalVector[0], normalVector[1], normalVector[2]);

src/mol-view/stage.ts

@@ -7,7 +7,7 @@
 import Viewer from './viewer'
 import { StateContext } from './state/context';
 import { Progress } from 'mol-task';
-import { MmcifUrlToModel, ModelToStructure, StructureToSpacefill, StructureToBallAndStick, StructureToDistanceRestraint, StructureToCartoon, StructureToBackbone } from './state/transform';
+import { MmcifUrlToModel, ModelToStructure, StructureToSpacefill, StructureToBallAndStick, StructureToDistanceRestraint, StructureToCartoon, StructureToBackbone, StructureCenter } from './state/transform';
 import { UrlEntity } from './state/entity';
 import { SpacefillProps } from 'mol-geo/representation/structure/spacefill';
 import { Context } from 'mol-app/context/context';
@@ -15,7 +15,7 @@ import { BallAndStickProps } from 'mol-geo/representation/structure/ball-and-sti
 import { CartoonProps } from 'mol-geo/representation/structure/cartoon';
 import { DistanceRestraintProps } from 'mol-geo/representation/structure/distance-restraint';
 import { BackboneProps } from 'mol-geo/representation/structure/backbone';
-import { Queries as Q, StructureProperties as SP, Query, Selection } from 'mol-model/structure';
+// import { Queries as Q, StructureProperties as SP, Query, Selection } from 'mol-model/structure';
 
 const spacefillProps: SpacefillProps = {
     doubleSided: true,
@@ -77,11 +77,13 @@ export class Stage {
         // this.loadPdbid('1hrv') // viral assembly
         // this.loadPdbid('1rb8') // virus
         // this.loadPdbid('1blu') // metal coordination
-        // this.loadPdbid('3pqr') // inter unit bonds
+        // this.loadPdbid('3pqr') // inter unit bonds, two polymer chains, ligands, water
         // this.loadPdbid('4v5a') // ribosome
         // this.loadPdbid('3j3q') // ...
         // this.loadPdbid('3sn6') // discontinuous chains
         // this.loadMmcifUrl(`../../examples/1cbs_full.bcif`)
+        // this.loadMmcifUrl(`../../examples/1cbs_updated.cif`)
+        // this.loadMmcifUrl(`../../examples/1crn.cif`)
 
         // this.loadMmcifUrl(`../../../test/pdb-dev/PDBDEV_00000001.cif`) // ok
         // this.loadMmcifUrl(`../../../test/pdb-dev/PDBDEV_00000002.cif`) // ok
@@ -101,24 +103,27 @@ export class Stage {
     async loadMmcifUrl (url: string) {
         const urlEntity = UrlEntity.ofUrl(this.ctx, url)
         const modelEntity = await MmcifUrlToModel.apply(this.ctx, urlEntity)
+        console.log(modelEntity.value)
         const structureEntity = await ModelToStructure.apply(this.ctx, modelEntity)
 
         StructureToBallAndStick.apply(this.ctx, structureEntity, { ...ballAndStickProps, visible: false })
         StructureToSpacefill.apply(this.ctx, structureEntity, { ...spacefillProps, visible: false })
         StructureToDistanceRestraint.apply(this.ctx, structureEntity, { ...distanceRestraintProps, visible: false })
-        // StructureToBackbone.apply(this.ctx, structureEntity, { ...backboneProps, visible: true })
+        StructureToBackbone.apply(this.ctx, structureEntity, { ...backboneProps, visible: true })
         StructureToCartoon.apply(this.ctx, structureEntity, { ...cartoonProps, visible: true })
+        StructureCenter.apply(this.ctx, structureEntity)
 
         this.globalContext.components.sequenceView.setState({ structure: structureEntity.value });
 
         // const structureEntity2 = await ModelToStructure.apply(this.ctx, modelEntity)
         // const q1 = Q.generators.atoms({
-        //     residueTest: l => SP.residue.label_seq_id(l) < 10
+        //     residueTest: l => SP.residue.label_seq_id(l) < 7
         // });
         // structureEntity2.value = Selection.unionStructure(await Query(q1)(structureEntity2.value).run());
-        // StructureToBackbone.apply(this.ctx, structureEntity2, { ...backboneProps, visible: true })
-        // StructureToCartoon.apply(this.ctx, structureEntity2, { ...cartoonProps, visible: true })
-        // StructureToBallAndStick.apply(this.ctx, structureEntity2, { ...ballAndStickProps, visible: true })
+        // await StructureToBackbone.apply(this.ctx, structureEntity2, { ...backboneProps, visible: true })
+        // await StructureToCartoon.apply(this.ctx, structureEntity2, { ...cartoonProps, visible: true })
+        // await StructureToBallAndStick.apply(this.ctx, structureEntity2, { ...ballAndStickProps, visible: false })
+        // StructureCenter.apply(this.ctx, structureEntity2)
     }
 
     loadPdbid (pdbid: string) {