1.2.2

- was ignored

package-lock.json

@@ -1,6 +1,6 @@
 {
   "name": "molstar",
-  "version": "1.2.0",
+  "version": "1.2.2",
   "lockfileVersion": 1,
   "requires": true,
   "dependencies": {

package.json

@@ -1,6 +1,6 @@
 {
   "name": "molstar",
-  "version": "1.2.0",
+  "version": "1.2.2",
   "description": "A comprehensive macromolecular library.",
   "homepage": "https://github.com/molstar/molstar#readme",
   "repository": {

src/apps/viewer/embedded.html

@@ -37,6 +37,7 @@
             });
             viewer.loadPdb('7bv2');
             viewer.loadEmdb('EMD-30210', { detail: 6 });
+            // viewer.loadAllModelsOrAssemblyFromUrl('https://cs.litemol.org/5ire/full', 'mmcif', false, { representationParams: { theme: { globalName: 'operator-name' } } })
         </script>
     </body>
 </html>

src/apps/viewer/index.ts

@@ -68,7 +68,7 @@ const DefaultViewerOptions = {
     layoutShowLog: true,
     layoutShowLeftPanel: true,
     disableAntialiasing: false,
-    pixelScale: void 0 as number | undefined,
+    pixelScale: 1,
     enableWboit: false,
 
     viewportShowExpand: PluginConfig.Viewport.ShowExpand.defaultValue,
@@ -115,16 +115,11 @@ export class Viewer {
             components: {
                 ...DefaultPluginSpec.components,
                 remoteState: o.layoutShowRemoteState ? 'default' : 'none',
-                viewport: {
-                    canvas3d: {
-                        multiSample: { mode: 'on', sampleLevel: 2  }
-                    }
-                }
             },
             config: [
                 [PluginConfig.General.DisableAntialiasing, o.disableAntialiasing],
                 [PluginConfig.General.PixelScale, o.pixelScale],
-                [PluginConfig.General.EnableWboit, true || o.enableWboit],
+                [PluginConfig.General.EnableWboit, o.enableWboit],
                 [PluginConfig.Viewport.ShowExpand, o.viewportShowExpand],
                 [PluginConfig.Viewport.ShowControls, o.viewportShowControls],
                 [PluginConfig.Viewport.ShowSettings, o.viewportShowSettings],
@@ -170,13 +165,13 @@ export class Viewer {
         }));
     }
 
-    async loadAllModelsOrAssemblyFromUrl(url: string, format: BuiltInTrajectoryFormat = 'mmcif', isBinary = false) {
+    async loadAllModelsOrAssemblyFromUrl(url: string, format: BuiltInTrajectoryFormat = 'mmcif', isBinary = false, options?: LoadStructureOptions) {
         const plugin = this.plugin;
 
         const data = await plugin.builders.data.download({ url, isBinary }, { state: { isGhost: true } });
         const trajectory = await plugin.builders.structure.parseTrajectory(data, format);
 
-        await this.plugin.builders.structure.hierarchy.applyPreset(trajectory, 'all-models', { useDefaultIfSingleModel: true });
+        await this.plugin.builders.structure.hierarchy.applyPreset(trajectory, 'all-models', { useDefaultIfSingleModel: true, representationPresetParams: options?.representationParams });
     }
 
     async loadStructureFromData(data: string | number[], format: BuiltInTrajectoryFormat, options?: { dataLabel?: string }) {

src/examples/alpha-orbitals/index.ts

@@ -4,9 +4,8 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
-import { Basis } from '../../extensions/alpha-orbitals/cubes';
-import { SphericalBasisOrder } from '../../extensions/alpha-orbitals/orbitals';
-import { CreateOrbitalRepresentation3D, CreateOrbitalVolume, StaticBasisAndOrbitals } from '../../extensions/alpha-orbitals/transforms';
+import { SphericalBasisOrder } from '../../extensions/alpha-orbitals/spherical-functions';
+import { BasisAndOrbitals, CreateOrbitalDensityVolume, CreateOrbitalRepresentation3D, CreateOrbitalVolume, StaticBasisAndOrbitals } from '../../extensions/alpha-orbitals/transforms';
 import { createPluginAsync, DefaultPluginSpec } from '../../mol-plugin';
 import { PluginStateObject } from '../../mol-plugin-state/objects';
 import { PluginConfig } from '../../mol-plugin/config';
@@ -20,23 +19,33 @@ import { DemoMoleculeSDF, DemoOrbitals } from './example-data';
 import { BehaviorSubject } from 'rxjs';
 import { debounceTime, skip } from 'rxjs/operators';
 import './index.html';
+import { Basis, AlphaOrbital } from '../../extensions/alpha-orbitals/data-model';
+import { canComputeAlphaOrbitalsOnGPU } from '../../extensions/alpha-orbitals/gpu/compute';
+import { PluginCommands } from '../../mol-plugin/commands';
 require('mol-plugin-ui/skin/light.scss');
 
 interface DemoInput {
     moleculeSdf: string,
     basis: Basis,
     order: SphericalBasisOrder,
-    orbitals: {
-        energy: number,
-        alpha: number[]
-    }[]
+    orbitals: AlphaOrbital[]
 }
 
 interface Params {
-    orbitalIndex: number,
+    show: { name: 'orbital', params: { index: number } } | { name: 'density', params: {} },
     isoValue: number,
-    gpuSurface: boolean,
-    cpuCompute: boolean
+    gpuSurface: boolean
+}
+
+type Selectors = {
+    type: 'orbital',
+    volume: StateObjectSelector<PluginStateObject.Volume.Data, typeof CreateOrbitalVolume>,
+    positive: StateObjectSelector<PluginStateObject.Volume.Representation3D, typeof CreateOrbitalRepresentation3D>
+    negative: StateObjectSelector<PluginStateObject.Volume.Representation3D, typeof CreateOrbitalRepresentation3D>
+} | {
+    type: 'density',
+    volume: StateObjectSelector<PluginStateObject.Volume.Data, typeof CreateOrbitalDensityVolume>,
+    positive: StateObjectSelector<PluginStateObject.Volume.Representation3D, typeof CreateOrbitalRepresentation3D>
 }
 
 export class AlphaOrbitalsExample {
@@ -62,6 +71,14 @@ export class AlphaOrbitalsExample {
 
         this.plugin.managers.interactivity.setProps({ granularity: 'element' });
 
+        if (!canComputeAlphaOrbitalsOnGPU(this.plugin.canvas3d?.webgl)) {
+            PluginCommands.Toast.Show(this.plugin, {
+                title: 'Error',
+                message: `Browser/device does not support required WebGL extension (OES_texture_float).`
+            });
+            return;
+        }
+
         this.load({
             moleculeSdf: DemoMoleculeSDF,
             ...DemoOrbitals
@@ -71,40 +88,84 @@ export class AlphaOrbitalsExample {
     }
 
     readonly params = new BehaviorSubject<ParamDefinition.For<Params>>({} as any);
-    readonly state = new BehaviorSubject<Params>({ orbitalIndex: 32, isoValue: 1, gpuSurface: false, cpuCompute: false });
+    readonly state = new BehaviorSubject<Params>({ show: { name: 'orbital', params: { index: 32 } }, isoValue: 1, gpuSurface: false });
+
+    private selectors?: Selectors = void 0;
+    private basis?: StateObjectSelector<BasisAndOrbitals> = void 0;
 
-    private volume?: StateObjectSelector<PluginStateObject.Volume.Data, typeof CreateOrbitalVolume>;
-    private positive?: StateObjectSelector<PluginStateObject.Volume.Representation3D, typeof CreateOrbitalRepresentation3D>;
-    private negative?: StateObjectSelector<PluginStateObject.Volume.Representation3D, typeof CreateOrbitalRepresentation3D>;
     private currentParams: Params = { ...this.state.value };
 
-    private async setIndex() {
-        if (!this.volume?.isOk) return;
+    private clearVolume() {
+        if (!this.selectors) return;
+        const v = this.selectors.volume;
+        this.selectors = void 0;
+        return this.plugin.build().delete(v).commit();
+    }
+
+    private async syncVolume() {
+        if (!this.basis?.isOk) return;
+
         const state = this.state.value;
-        await this.plugin.build().to(this.volume).update(CreateOrbitalVolume, () => ({ index: state.orbitalIndex, cpuCompute: state.cpuCompute })).commit();
+
+        if (state.show.name !== this.selectors?.type) {
+            await this.clearVolume();
+        }
+
+        const update = this.plugin.build();
+        if (state.show.name === 'orbital') {
+            if (!this.selectors) {
+                const volume = update
+                    .to(this.basis)
+                    .apply(CreateOrbitalVolume, { index: state.show.params.index });
+
+                const positive = volume.apply(CreateOrbitalRepresentation3D, this.volumeParams('positive', ColorNames.blue)).selector;
+                const negative = volume.apply(CreateOrbitalRepresentation3D, this.volumeParams('negative', ColorNames.red)).selector;
+
+                this.selectors = { type: 'orbital', volume: volume.selector, positive, negative };
+            } else {
+                const index = state.show.params.index;
+                update.to(this.selectors.volume).update(CreateOrbitalVolume, () => ({ index }));
+            }
+        } else {
+            if (!this.selectors) {
+                const volume = update
+                    .to(this.basis)
+                    .apply(CreateOrbitalDensityVolume);
+                const positive = volume.apply(CreateOrbitalRepresentation3D, this.volumeParams('positive', ColorNames.blue)).selector;
+                this.selectors = { type: 'density', volume: volume.selector, positive };
+            }
+        }
+
+        await update.commit();
+
         if (this.currentParams.gpuSurface !== this.state.value.gpuSurface) {
             await this.setIsovalue();
         }
+
         this.currentParams = this.state.value;
     }
 
     private setIsovalue() {
+        if (!this.selectors) return;
+
         this.currentParams = this.state.value;
         const update = this.plugin.build();
-        update.to(this.positive!).update(this.volumeParams('positive', ColorNames.blue));
-        update.to(this.negative!).update(this.volumeParams('negative', ColorNames.red));
+        update.to(this.selectors.positive).update(this.volumeParams('positive', ColorNames.blue));
+        if (this.selectors?.type === 'orbital') {
+            update.to(this.selectors.negative).update(this.volumeParams('negative', ColorNames.red));
+        }
         return update.commit();
     }
 
     private volumeParams(kind: 'positive' | 'negative', color: Color): StateTransformer.Params<typeof CreateOrbitalRepresentation3D> {
         return {
-            alpha: 1.0,
+            alpha: 0.85,
             color,
             directVolume: this.state.value.gpuSurface,
             kind,
             relativeIsovalue: this.state.value.isoValue,
             pickable: false,
-            xrayShaded: false
+            xrayShaded: true
         };
     }
 
@@ -119,36 +180,28 @@ export class AlphaOrbitalsExample {
         const all = await this.plugin.builders.structure.tryCreateComponentStatic(structure, 'all');
         if (all) await this.plugin.builders.structure.representation.addRepresentation(all, { type: 'ball-and-stick', color: 'element-symbol', colorParams: { carbonColor: { name: 'element-symbol', params: {} } } });
 
-        const state = this.state.value;
 
-        this.volume = await this.plugin.build().toRoot()
+        this.basis = await this.plugin.build().toRoot()
             .apply(StaticBasisAndOrbitals, { basis: input.basis, order: input.order, orbitals: input.orbitals })
-            .apply(CreateOrbitalVolume, { index: state.orbitalIndex, forceCpuCompute: this.currentParams.cpuCompute })
             .commit();
 
-        if (!this.volume.isOk) {
-            this.volume = void 0;
-            return;
-        }
-
-        const repr = this.plugin.build().to(this.volume);
-
-        this.positive = repr.apply(CreateOrbitalRepresentation3D, this.volumeParams('positive', ColorNames.blue)).selector;
-        this.negative = repr.apply(CreateOrbitalRepresentation3D, this.volumeParams('negative', ColorNames.red)).selector;
-
-        await repr.commit();
+        await this.syncVolume();
 
         this.params.next({
-            orbitalIndex: ParamDefinition.Numeric(this.currentParams.orbitalIndex, { min: 0, max: input.orbitals.length - 1 }, { immediateUpdate: true, isEssential: true }),
+            show: ParamDefinition.MappedStatic('orbital', {
+                'orbital': ParamDefinition.Group({
+                    index: ParamDefinition.Numeric(32, { min: 0, max: input.orbitals.length - 1 }, { immediateUpdate: true, isEssential: true }),
+                }),
+                'density': ParamDefinition.EmptyGroup()
+            }, { cycle: true }),
             isoValue: ParamDefinition.Numeric(this.currentParams.isoValue, { min: 0.5, max: 3, step: 0.1 }, { immediateUpdate: true, isEssential: false }),
-            gpuSurface: ParamDefinition.Boolean(this.currentParams.gpuSurface),
-            cpuCompute: ParamDefinition.Boolean(this.currentParams.cpuCompute, { label: 'CPU Compute' })
+            gpuSurface: ParamDefinition.Boolean(this.currentParams.gpuSurface, { isHidden: true })
         });
 
         this.state.pipe(skip(1), debounceTime(1000 / 24)).subscribe(async params => {
-            if (params.orbitalIndex !== this.currentParams.orbitalIndex
-                || params.cpuCompute !== this.currentParams.cpuCompute) {
-                this.setIndex();
+            if (params.show.name !== this.currentParams.show.name
+                || (params.show.name === 'orbital' && this.currentParams.show.name === 'orbital' && params.show.params.index !== this.currentParams.show.params.index)) {
+                this.syncVolume();
             } else if (params.isoValue !== this.currentParams.isoValue || params.gpuSurface !== this.currentParams.gpuSurface) {
                 this.setIsovalue();
             }

src/examples/alpha-orbitals/test-water.ts

@@ -1,190 +0,0 @@
-import { CollocationParams } from '../../extensions/alpha-orbitals/collocation';
-import { Basis } from '../../extensions/alpha-orbitals/cubes';
-import { Box3D } from '../../mol-math/geometry';
-import { Vec3 } from '../../mol-math/linear-algebra';
-
-const _testBasis: Basis = {
-    'atoms': [
-        {
-            'center': [
-                0.025886090588624934,
-                0.019164790004065606,
-                -0.013539970104105408
-            ] as Vec3,
-            'shells': [
-                {
-                    'angularMomentum': [0],
-                    'coefficients': [
-                        [
-                            -0.004151277818987536,
-                            -0.02067024147993795,
-                            -0.05150303336984537,
-                            0.33462711739899537,
-                            0.5621061300983125,
-                            0.17129946969948573
-                        ]
-                    ],
-                    'exponents': [
-                        152.28769660788095,
-                        27.928015215973073,
-                        7.848374792384515,
-                        1.1223350202705642,
-                        0.5093846587907856,
-                        0.24292266532510307
-                    ]
-                },
-                {
-                    'angularMomentum': [1],
-                    'coefficients': [
-                        [
-                            0.007924233646294425,
-                            0.051441048251911314,
-                            0.18984000600705359,
-                            0.4049863191150474,
-                            0.40123628611490797,
-                            0.1051855189039082
-                        ]
-                    ],
-                    'exponents': [
-                        27.203421487167727,
-                        7.09409912597673,
-                        2.5383362605345954,
-                        1.0610730767843852,
-                        0.4851948916410433,
-                        0.22938302550642545
-                    ]
-                }
-            ]
-        },
-        {
-            'center': [
-                0.5082729578468134,
-                1.6880351220025265,
-                0.4963443067810461
-            ] as Vec3,
-            'shells': [
-                {
-                    'angularMomentum': [0],
-                    'coefficients': [
-                        [
-                            0.009163596280542963,
-                            0.04936149294292479,
-                            0.16853830490998634,
-                            0.37056279972195677,
-                            0.4164915298246781,
-                            0.13033408410772263
-                        ]
-                    ],
-                    'exponents': [
-                        33.710073211949485,
-                        6.180705022740464,
-                        1.7291385346152253,
-                        0.5940057549921978,
-                        0.2306698170449518,
-                        0.09500256906284119
-                    ]
-                },
-                {
-                    'angularMomentum': [0],
-                    'coefficients': [
-                        [
-                            -0.32279868167000036,
-                            3.209629817295221,
-                            2.4672629224617935,
-                            -0.048487066612842224,
-                            -0.2611850111200143,
-                            -0.8917817597810863,
-                            -1.9607480081275706,
-                            -2.203769342520311,
-                            -0.6896328935259993
-                        ]
-                    ],
-                    'exponents': [
-                        10.256286070314905,
-                        0.6227965325875392,
-                        0.2391007667853915,
-                        33.710073211949485,
-                        6.180705022740464,
-                        1.7291385346152253,
-                        0.5940057549921978,
-                        0.2306698170449518,
-                        0.09500256906284119
-                    ]
-                }
-            ]
-        },
-        {
-            'center': [
-                1.1367367844436005,
-                -0.47018519422670163,
-                -1.356802622574504
-            ] as Vec3,
-            'shells': [
-                {
-                    'angularMomentum': [0],
-                    'coefficients': [
-                        [
-                            0.009163596280542963,
-                            0.04936149294292479,
-                            0.16853830490998634,
-                            0.37056279972195677,
-                            0.4164915298246781,
-                            0.13033408410772263
-                        ]
-                    ],
-                    'exponents': [
-                        33.710073211949485,
-                        6.180705022740464,
-                        1.7291385346152253,
-                        0.5940057549921978,
-                        0.2306698170449518,
-                        0.09500256906284119
-                    ]
-                },
-                {
-                    'angularMomentum': [0],
-                    'coefficients': [
-                        [
-                            -0.32279868167000036,
-                            3.209629817295221,
-                            2.4672629224617935,
-                            -0.048487066612842224,
-                            -0.2611850111200143,
-                            -0.8917817597810863,
-                            -1.9607480081275706,
-                            -2.203769342520311,
-                            -0.6896328935259993
-                        ]
-                    ],
-                    'exponents': [
-                        10.256286070314905,
-                        0.6227965325875392,
-                        0.2391007667853915,
-                        33.710073211949485,
-                        6.180705022740464,
-                        1.7291385346152253,
-                        0.5940057549921978,
-                        0.2306698170449518,
-                        0.09500256906284119
-                    ]
-                }
-            ]
-        }
-    ]
-};
-
-const grid = {
-    box: Box3D.create(Vec3.create(-1, -1, -1), Vec3.create(1, 1, 1)),
-    delta: Vec3.create(2, 2, 2),
-    dimensions: Vec3.create(2, 2, 2),
-    npoints: 8,
-    size: Vec3.create(2, 2, 2)
-};
-
-export const TestWaterParams: CollocationParams = {
-    grid,
-    alphaOrbitals: [-2.2623991420609075e-16, 0.6360205395000592, 0.6672122399886391, -0.3876927909355508, -1.6780131293332933e-16, 2.844782862661151e-16, 4.977960694176068e-19, -2.3945919908996803e-16],
-    basis: _testBasis,
-    cutoffThreshold: 0,
-    sphericalOrder: 'cca-reverse'
-};

src/extensions/alpha-orbitals/_spec/collocation.spec.ts

@@ -8,11 +8,18 @@ import { Box3D } from '../../../mol-math/geometry';
 import { Vec3 } from '../../../mol-math/linear-algebra';
 import { RuntimeContext } from '../../../mol-task';
 import { sphericalCollocation } from '../collocation';
-import { Basis } from '../cubes';
+import { Basis, CubeGridInfo } from '../data-model';
 
 describe('alpha-orbitals-cubes', () => {
     it('water', async () => {
-        const grid = {
+        const grid: CubeGridInfo = {
+            params: {
+                basis: _testBasis,
+                cutoffThreshold: 0,
+                sphericalOrder: 'cca-reverse',
+                boxExpand: 0,
+                gridSpacing: []
+            },
             box: Box3D.create(Vec3.create(-1, -1, -1), Vec3.create(1, 1, 1)),
             delta: Vec3.create(2, 2, 2),
             dimensions: Vec3.create(2, 2, 2),
@@ -20,12 +27,10 @@ describe('alpha-orbitals-cubes', () => {
             size: Vec3.create(2, 2, 2)
         };
 
-        const matrix = await sphericalCollocation({
-            grid,
-            alphaOrbitals: [-2.2623991420609075e-16, 0.6360205395000592, 0.6672122399886391, -0.3876927909355508, -1.6780131293332933e-16, 2.844782862661151e-16, 4.977960694176068e-19, -2.3945919908996803e-16],
-            basis: _testBasis,
-            cutoffThreshold: 0,
-            sphericalOrder: 'cca-reverse'
+        const matrix = await sphericalCollocation(grid, {
+            energy: 0,
+            occupancy: 0,
+            alpha: [-2.2623991420609075e-16, 0.6360205395000592, 0.6672122399886391, -0.3876927909355508, -1.6780131293332933e-16, 2.844782862661151e-16, 4.977960694176068e-19, -2.3945919908996803e-16]
         }, RuntimeContext.Synchronous);
 
         const expected = [-0.1451730622877498, 0.06479453956039086, -0.2777738736440713, -0.057116584776260436, 0.05929916178822645, 0.2742903371231049, -0.07221698722165386, 0.15389180241391376];

src/extensions/alpha-orbitals/collocation.ts

@@ -9,31 +9,15 @@
 import { Vec3 } from '../../mol-math/linear-algebra';
 import { RuntimeContext } from '../../mol-task';
 import { arrayMin } from '../../mol-util/array';
-import { Basis, CubeGridInfo } from './cubes';
-import {
-    normalizeBasicOrder,
-
-    SphericalBasisOrder, SphericalFunctions
-} from './orbitals';
-
-export interface CollocationParams {
-    grid: CubeGridInfo;
-    basis: Basis;
-    sphericalOrder: SphericalBasisOrder;
-    cutoffThreshold: number;
-    alphaOrbitals: number[];
-}
+import { AlphaOrbital, CubeGridInfo } from './data-model';
+import { normalizeBasicOrder, SphericalFunctions } from './spherical-functions';
 
 export async function sphericalCollocation(
-    {
-        grid,
-        basis,
-        sphericalOrder,
-        cutoffThreshold,
-        alphaOrbitals,
-    }: CollocationParams,
+    grid: CubeGridInfo,
+    orbital: AlphaOrbital,
     taskCtx: RuntimeContext
 ) {
+    const { basis, sphericalOrder, cutoffThreshold } = grid.params;
     let baseCount = 0;
 
     for (const atom of basis.atoms) {
@@ -58,7 +42,7 @@ export async function sphericalCollocation(
             for (const L of shell.angularMomentum) {
                 const alpha = normalizeBasicOrder(
                     L,
-                    alphaOrbitals.slice(baseIndex, baseIndex + 2 * L + 1),
+                    orbital.alpha.slice(baseIndex, baseIndex + 2 * L + 1),
                     sphericalOrder
                 );
                 baseIndex += 2 * L + 1;

src/extensions/alpha-orbitals/cubes.ts

@@ -1,260 +0,0 @@
-/**
- * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * Inspired by https://github.com/dgasmith/gau2grid.
- *
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-
-import { sortArray } from '../../mol-data/util';
-import { WebGLContext } from '../../mol-gl/webgl/context';
-import { Box3D } from '../../mol-math/geometry';
-import { Mat4, Tensor, Vec3 } from '../../mol-math/linear-algebra';
-import { Grid } from '../../mol-model/volume';
-import { Task } from '../../mol-task';
-import { arrayMax, arrayMin, arrayRms } from '../../mol-util/array';
-import { CollocationParams, sphericalCollocation } from './collocation';
-import { canComputeAlphaOrbitalsOnGPU, gpuComputeAlphaOrbitalsGridValues } from './gpu/compute';
-import { SphericalBasisOrder } from './orbitals';
-
-export interface CubeGridInfo {
-    dimensions: Vec3;
-    box: Box3D;
-    size: Vec3;
-    npoints: number;
-    delta: Vec3;
-}
-
-export interface CubeGrid {
-    grid: Grid;
-    isovalues?: { negative?: number; positive?: number };
-}
-
-export interface Basis {
-    atoms: {
-        // in Bohr units!
-        center: Vec3;
-        shells: SphericalElectronShell[];
-    }[];
-}
-
-// Note: generally contracted gaussians are currently not supported.
-export interface SphericalElectronShell {
-    exponents: number[];
-    angularMomentum: number[];
-    // number[] for each angular momentum
-    coefficients: number[][];
-}
-
-export interface SphericalCollocationParams {
-    basis: Basis;
-    /**
-     * for each electron shell compute a cutoff radius as
-     *
-     *    const cutoffRadius = Math.sqrt(-Math.log(cutoffThreshold) / arrayMin(exponents));
-     *
-     */
-    cutoffThreshold: number;
-    sphericalOrder: SphericalBasisOrder;
-    boxExpand: number;
-    gridSpacing: number | [atomCountThreshold: number, spacing: number][];
-    alphaOrbitals: number[];
-    doNotComputeIsovalues?: boolean
-}
-
-export function createSphericalCollocationGrid(
-    params: SphericalCollocationParams, webgl?: WebGLContext
-): Task<CubeGrid> {
-    return Task.create('Spherical Collocation Grid', async (ctx) => {
-        const centers = params.basis.atoms.map(a => a.center);
-
-        const cParams: CollocationParams = {
-            grid: initBox(centers, params.gridSpacing, params.boxExpand),
-            basis: params.basis,
-            alphaOrbitals: params.alphaOrbitals,
-            cutoffThreshold: params.cutoffThreshold,
-            sphericalOrder: params.sphericalOrder
-        };
-
-        let matrix: Float32Array;
-        if (canComputeAlphaOrbitalsOnGPU(webgl)) {
-            // console.time('gpu');
-            matrix = gpuComputeAlphaOrbitalsGridValues(webgl!, cParams);
-            // console.timeEnd('gpu');
-        } else {
-            // console.time('cpu');
-            matrix = await sphericalCollocation(cParams, ctx);
-            // console.timeEnd('cpu');
-        }
-
-        return createCubeGrid(cParams.grid, matrix, [0, 1, 2], !params.doNotComputeIsovalues);
-    });
-}
-
-const BohrToAngstromFactor = 0.529177210859;
-
-function createCubeGrid(gridInfo: CubeGridInfo, values: Float32Array, axisOrder: number[], computeIsovalues: boolean) {
-    const boxSize = Box3D.size(Vec3(), gridInfo.box);
-    const boxOrigin = Vec3.clone(gridInfo.box.min);
-
-    Vec3.scale(boxSize, boxSize, BohrToAngstromFactor);
-    Vec3.scale(boxOrigin, boxOrigin, BohrToAngstromFactor);
-
-    const scale = Mat4.fromScaling(
-        Mat4(),
-        Vec3.div(
-            Vec3(),
-            boxSize,
-            Vec3.sub(Vec3(), gridInfo.dimensions, Vec3.create(1, 1, 1))
-        )
-    );
-    const translate = Mat4.fromTranslation(Mat4(), boxOrigin);
-    const matrix = Mat4.mul(Mat4(), translate, scale);
-
-    const grid: Grid = {
-        transform: { kind: 'matrix', matrix },
-        cells: Tensor.create(
-            Tensor.Space(gridInfo.dimensions, axisOrder, Float32Array),
-            (values as any) as Tensor.Data
-        ),
-        stats: {
-            min: arrayMin(values),
-            max: arrayMax(values),
-            mean: arrayMax(values),
-            sigma: arrayRms(values),
-        },
-    };
-
-    // TODO: when using GPU rendering, the cumulative sum can be computed
-    //       along the ray on the fly?
-    let isovalues: { negative?: number, positive?: number } | undefined;
-
-    if (computeIsovalues) {
-        isovalues = computeIsocontourValues(values, 0.85);
-    }
-
-    return { grid, isovalues };
-}
-
-function initBox(
-    geometry: Vec3[],
-    spacing: SphericalCollocationParams['gridSpacing'],
-    expand: number
-): CubeGridInfo {
-    const count = geometry.length;
-    const box = Box3D.expand(
-        Box3D(),
-        Box3D.fromVec3Array(Box3D(), geometry),
-        Vec3.create(expand, expand, expand)
-    );
-    const size = Box3D.size(Vec3(), box);
-
-    const spacingThresholds =
-        typeof spacing === 'number' ? [[0, spacing]] : [...spacing];
-    spacingThresholds.sort((a, b) => b[0] - a[0]);
-
-    let s = 0.4;
-    for (let i = 0; i <= spacingThresholds.length; i++) {
-        s = spacingThresholds[i][1];
-        if (spacingThresholds[i][0] <= count) break;
-    }
-
-    const dimensions = Vec3.ceil(Vec3(), Vec3.scale(Vec3(), size, 1 / s));
-
-    return {
-        box,
-        dimensions,
-        size,
-        npoints: dimensions[0] * dimensions[1] * dimensions[2],
-        delta: Vec3.div(Vec3(), size, Vec3.subScalar(Vec3(), dimensions, 1)),
-    };
-}
-
-export function computeIsocontourValues(
-    input: Float32Array,
-    cumulativeThreshold: number
-) {
-    let weightSum = 0;
-    for (let i = 0, _i = input.length; i < _i; i++) {
-        const v = input[i];
-        const w = v * v;
-        weightSum += w;
-    }
-    const avgWeight = weightSum / input.length;
-    let minWeight = 3 * avgWeight;
-
-    // do not try to identify isovalues for degenerate data
-    // e.g. all values are almost same
-    if (Math.abs(avgWeight - input[0] * input[0]) < 1e-5) {
-        return { negative: void 0, positive: void 0 };
-    }
-
-    let size = 0;
-    while (true) {
-        let csum = 0;
-        size = 0;
-        for (let i = 0, _i = input.length; i < _i; i++) {
-            const v = input[i];
-            const w = v * v;
-            if (w >= minWeight) {
-                csum += w;
-                size++;
-            }
-        }
-
-        if (csum / weightSum > cumulativeThreshold) {
-            break;
-        }
-
-        minWeight -= avgWeight;
-    }
-
-    const values = new Float32Array(size);
-    const weights = new Float32Array(size);
-    const indices = new Int32Array(size);
-
-    let o = 0;
-    for (let i = 0, _i = input.length; i < _i; i++) {
-        const v = input[i];
-        const w = v * v;
-        if (w >= minWeight) {
-            values[o] = v;
-            weights[o] = w;
-            indices[o] = o;
-            o++;
-        }
-    }
-
-    sortArray(
-        indices,
-        (indices, i, j) => weights[indices[j]] - weights[indices[i]]
-    );
-
-    let cweight = 0,
-        cutoffIndex = 0;
-    for (let i = 0; i < size; i++) {
-        cweight += weights[indices[i]];
-
-        if (cweight / weightSum >= cumulativeThreshold) {
-            cutoffIndex = i;
-            break;
-        }
-    }
-
-    let positive = Number.POSITIVE_INFINITY,
-        negative = Number.NEGATIVE_INFINITY;
-
-    for (let i = 0; i < cutoffIndex; i++) {
-        const v = values[indices[i]];
-        if (v > 0) {
-            if (v < positive) positive = v;
-        } else if (v < 0) {
-            if (v > negative) negative = v;
-        }
-    }
-
-    return {
-        negative: negative !== Number.NEGATIVE_INFINITY ? negative : void 0,
-        positive: positive !== Number.POSITIVE_INFINITY ? positive : void 0,
-    };
-}

src/extensions/alpha-orbitals/data-model.ts

@@ -0,0 +1,131 @@
+/**
+ * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Mat4, Tensor, Vec3 } from '../../mol-math/linear-algebra';
+import { Grid } from '../../mol-model/volume';
+import { SphericalBasisOrder } from './spherical-functions';
+import { Box3D } from '../../mol-math/geometry';
+import { arrayMin, arrayMax, arrayRms } from '../../mol-util/array';
+
+// Note: generally contracted gaussians are currently not supported.
+export interface SphericalElectronShell {
+    exponents: number[];
+    angularMomentum: number[];
+    // number[] for each angular momentum
+    coefficients: number[][];
+}
+
+export interface Basis {
+    atoms: {
+        // in Bohr units!
+        center: Vec3;
+        shells: SphericalElectronShell[];
+    }[];
+}
+
+export interface AlphaOrbital {
+    energy: number;
+    occupancy: number;
+    alpha: number[];
+}
+
+export interface CubeGridComputationParams {
+    basis: Basis;
+    /**
+     * for each electron shell compute a cutoff radius as
+     *    const cutoffRadius = Math.sqrt(-Math.log(cutoffThreshold) / arrayMin(exponents));
+     */
+    cutoffThreshold: number;
+    sphericalOrder: SphericalBasisOrder;
+    boxExpand: number;
+    gridSpacing: number | [atomCountThreshold: number, spacing: number][];
+    doNotComputeIsovalues?: boolean;
+}
+
+export interface CubeGridInfo {
+    params: CubeGridComputationParams;
+    dimensions: Vec3;
+    box: Box3D;
+    size: Vec3;
+    npoints: number;
+    delta: Vec3;
+}
+
+export interface CubeGrid {
+    grid: Grid;
+    isovalues?: { negative?: number; positive?: number };
+}
+
+export function initCubeGrid(params: CubeGridComputationParams): CubeGridInfo {
+    const geometry = params.basis.atoms.map(a => a.center);
+    const { gridSpacing: spacing, boxExpand: expand } = params;
+
+    const count = geometry.length;
+    const box = Box3D.expand(
+        Box3D(),
+        Box3D.fromVec3Array(Box3D(), geometry),
+        Vec3.create(expand, expand, expand)
+    );
+    const size = Box3D.size(Vec3(), box);
+
+    const spacingThresholds =
+        typeof spacing === 'number' ? [[0, spacing]] : [...spacing];
+    spacingThresholds.sort((a, b) => b[0] - a[0]);
+
+    let s = 0.4;
+    for (let i = 0; i <= spacingThresholds.length; i++) {
+        s = spacingThresholds[i][1];
+        if (spacingThresholds[i][0] <= count) break;
+    }
+
+    const dimensions = Vec3.ceil(Vec3(), Vec3.scale(Vec3(), size, 1 / s));
+
+    return {
+        params,
+        box,
+        dimensions,
+        size,
+        npoints: dimensions[0] * dimensions[1] * dimensions[2],
+        delta: Vec3.div(Vec3(), size, Vec3.subScalar(Vec3(), dimensions, 1)),
+    };
+}
+
+const BohrToAngstromFactor = 0.529177210859;
+
+export function createGrid(gridInfo: CubeGridInfo, values: Float32Array, axisOrder: number[]) {
+    const boxSize = Box3D.size(Vec3(), gridInfo.box);
+    const boxOrigin = Vec3.clone(gridInfo.box.min);
+
+    Vec3.scale(boxSize, boxSize, BohrToAngstromFactor);
+    Vec3.scale(boxOrigin, boxOrigin, BohrToAngstromFactor);
+
+    const scale = Mat4.fromScaling(
+        Mat4(),
+        Vec3.div(
+            Vec3(),
+            boxSize,
+            Vec3.sub(Vec3(), gridInfo.dimensions, Vec3.create(1, 1, 1))
+        )
+    );
+    const translate = Mat4.fromTranslation(Mat4(), boxOrigin);
+    const matrix = Mat4.mul(Mat4(), translate, scale);
+
+    const grid: Grid = {
+        transform: { kind: 'matrix', matrix },
+        cells: Tensor.create(
+            Tensor.Space(gridInfo.dimensions, axisOrder, Float32Array),
+            (values as any) as Tensor.Data
+        ),
+        stats: {
+            min: arrayMin(values),
+            max: arrayMax(values),
+            mean: arrayMax(values),
+            sigma: arrayRms(values),
+        },
+    };
+
+    return grid;
+}

src/extensions/alpha-orbitals/density.ts

@@ -0,0 +1,123 @@
+/**
+ * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { sortArray } from '../../mol-data/util';
+import { WebGLContext } from '../../mol-gl/webgl/context';
+import { Task } from '../../mol-task';
+import { AlphaOrbital, createGrid, CubeGrid, CubeGridComputationParams, initCubeGrid } from './data-model';
+import { canComputeAlphaOrbitalsOnGPU, gpuComputeAlphaOrbitalsDensityGridValues } from './gpu/compute';
+
+export function createSphericalCollocationDensityGrid(
+    params: CubeGridComputationParams, orbitals: AlphaOrbital[], webgl?: WebGLContext
+): Task<CubeGrid> {
+    return Task.create('Spherical Collocation Grid', async (ctx) => {
+        const cubeGrid = initCubeGrid(params);
+
+        let matrix: Float32Array;
+        if (canComputeAlphaOrbitalsOnGPU(webgl)) {
+            // console.time('gpu');
+            matrix = await gpuComputeAlphaOrbitalsDensityGridValues(webgl!, cubeGrid, orbitals, ctx);
+            // console.timeEnd('gpu');
+        } else {
+            throw new Error('Missing OES_texture_float WebGL extension.');
+        }
+
+        const grid = createGrid(cubeGrid, matrix, [0, 1, 2]);
+        let isovalues: { negative?: number, positive?: number } | undefined;
+
+        if (!params.doNotComputeIsovalues) {
+            isovalues = computeDensityIsocontourValues(matrix, 0.85);
+        }
+
+        return { grid, isovalues };
+    });
+}
+
+export function computeDensityIsocontourValues(input: Float32Array, cumulativeThreshold: number) {
+    let weightSum = 0;
+    for (let i = 0, _i = input.length; i < _i; i++) {
+        const v = input[i];
+        const w = Math.abs(v);
+        weightSum += w;
+    }
+    const avgWeight = weightSum / input.length;
+    let minWeight = 3 * avgWeight;
+
+    // do not try to identify isovalues for degenerate data
+    // e.g. all values are almost same
+    if (Math.abs(avgWeight - input[0] * input[0]) < 1e-5) {
+        return { negative: void 0, positive: void 0 };
+    }
+
+    let size = 0;
+    while (true) {
+        let csum = 0;
+        size = 0;
+        for (let i = 0, _i = input.length; i < _i; i++) {
+            const v = input[i];
+            const w = Math.abs(v);
+            if (w >= minWeight) {
+                csum += w;
+                size++;
+            }
+        }
+
+        if (csum / weightSum > cumulativeThreshold) {
+            break;
+        }
+
+        minWeight -= avgWeight;
+    }
+
+    const values = new Float32Array(size);
+    const weights = new Float32Array(size);
+    const indices = new Int32Array(size);
+
+    let o = 0;
+    for (let i = 0, _i = input.length; i < _i; i++) {
+        const v = input[i];
+        const w = Math.abs(v);
+        if (w >= minWeight) {
+            values[o] = v;
+            weights[o] = w;
+            indices[o] = o;
+            o++;
+        }
+    }
+
+    sortArray(
+        indices,
+        (indices, i, j) => weights[indices[j]] - weights[indices[i]]
+    );
+
+    let cweight = 0,
+        cutoffIndex = 0;
+    for (let i = 0; i < size; i++) {
+        cweight += weights[indices[i]];
+
+        if (cweight / weightSum >= cumulativeThreshold) {
+            cutoffIndex = i;
+            break;
+        }
+    }
+
+    let positive = Number.POSITIVE_INFINITY,
+        negative = Number.NEGATIVE_INFINITY;
+
+    for (let i = 0; i < cutoffIndex; i++) {
+        const v = values[indices[i]];
+        if (v > 0) {
+            if (v < positive) positive = v;
+        } else if (v < 0) {
+            if (v > negative) negative = v;
+        }
+    }
+
+    return {
+        negative: negative !== Number.NEGATIVE_INFINITY ? negative : void 0,
+        positive: positive !== Number.POSITIVE_INFINITY ? positive : void 0,
+    };
+}

src/extensions/alpha-orbitals/gpu/compute.ts

@@ -11,11 +11,12 @@ import { ShaderCode } from '../../../mol-gl/shader-code';
 import quad_vert from '../../../mol-gl/shader/quad.vert';
 import { WebGLContext } from '../../../mol-gl/webgl/context';
 import { createComputeRenderItem } from '../../../mol-gl/webgl/render-item';
+import { RuntimeContext } from '../../../mol-task';
 import { ValueCell } from '../../../mol-util';
 import { arrayMin } from '../../../mol-util/array';
 import { isLittleEndian } from '../../../mol-util/is-little-endian';
-import { CollocationParams } from '../collocation';
-import { normalizeBasicOrder } from '../orbitals';
+import { AlphaOrbital, Basis, CubeGridInfo } from '../data-model';
+import { normalizeBasicOrder, SphericalBasisOrder } from '../spherical-functions';
 import shader_frag from './shader.frag';
 
 const AlphaOrbitalsSchema = {
@@ -32,18 +33,38 @@ const AlphaOrbitalsSchema = {
     uNAlpha: UniformSpec('i'),
     uNCoeff: UniformSpec('i'),
     uMaxCoeffs: UniformSpec('i'),
-    uLittleEndian: UniformSpec('b')
+    uLittleEndian: UniformSpec('b'),
+    uDensity: UniformSpec('b'),
+    uOccupancy: UniformSpec('f'),
+    tCumulativeSum: TextureSpec('texture', 'rgba', 'ubyte', 'nearest')
 };
+type AlphaOrbitalsSchema = Values<typeof AlphaOrbitalsSchema>
 const AlphaOrbitalsName = 'alpha-orbitals';
+const AlphaOrbitalsTex0 = 'alpha-orbitals-0';
+const AlphaOrbitalsTex1 = 'alpha-orbitals-1';
 const AlphaOrbitalsShaderCode = ShaderCode(AlphaOrbitalsName, quad_vert, shader_frag);
-type AlphaOrbitalsRenderable = ComputeRenderable<Values<typeof AlphaOrbitalsSchema>>
+type AlphaOrbitalsRenderable = ComputeRenderable<AlphaOrbitalsSchema>
+
+function getNormalizedAlpha(basis: Basis, alphaOrbitals: number[], sphericalOrder: SphericalBasisOrder) {
+    const alpha = new Float32Array(alphaOrbitals.length);
+
+    let aO = 0;
+    for (const atom of basis.atoms) {
+        for (const shell of atom.shells) {
+            for (const L of shell.angularMomentum) {
+                const a0 = normalizeBasicOrder(L, alphaOrbitals.slice(aO, aO + 2 * L + 1), sphericalOrder);
+                for (let i = 0; i < a0.length; i++) alpha[aO + i] = a0[i];
+                aO += 2 * L + 1;
+            }
+        }
+    }
+
+    return alpha;
+}
+
+function createTextureData(grid: CubeGridInfo, orbital: AlphaOrbital) {
+    const { basis, sphericalOrder, cutoffThreshold } = grid.params;
 
-function createTextureData({
-    basis,
-    sphericalOrder,
-    alphaOrbitals,
-    cutoffThreshold
-}: CollocationParams) {
     let centerCount = 0;
     let baseCount = 0;
     let coeffCount = 0;
@@ -75,7 +96,7 @@ function createTextureData({
 
             let amIndex = 0;
             for (const L of shell.angularMomentum) {
-                const a0 = normalizeBasicOrder(L, alphaOrbitals.slice(aO, aO + 2 * L + 1), sphericalOrder);
+                const a0 = normalizeBasicOrder(L, orbital.alpha.slice(aO, aO + 2 * L + 1), sphericalOrder);
 
                 const cutoffRadius = cutoffThreshold > 0
                     ? Math.sqrt(-Math.log(cutoffThreshold) / arrayMin(shell.exponents))
@@ -113,17 +134,27 @@ function createTextureData({
     return { nCenters: centerCount, nAlpha: baseCount, nCoeff: coeffCount, maxCoeffs, centers, info, alpha, coeff };
 }
 
-function createAlphaOrbitalsRenderable(ctx: WebGLContext, params: CollocationParams): AlphaOrbitalsRenderable {
-    const data = createTextureData(params);
+function createAlphaOrbitalsRenderable(ctx: WebGLContext, grid: CubeGridInfo, orbital: AlphaOrbital): AlphaOrbitalsRenderable {
+    const data = createTextureData(grid, orbital);
 
-    const [nx, ny, nz] = params.grid.dimensions;
+    const [nx, ny, nz] = grid.dimensions;
     const width = Math.ceil(Math.sqrt(nx * ny * nz));
 
-    const values: Values<typeof AlphaOrbitalsSchema> = {
+    if (!ctx.namedFramebuffers[AlphaOrbitalsName]) {
+        ctx.namedFramebuffers[AlphaOrbitalsName] = ctx.resources.framebuffer();
+    }
+    if (!ctx.namedTextures[AlphaOrbitalsTex0]) {
+        ctx.namedTextures[AlphaOrbitalsTex0] = ctx.resources.texture('image-uint8', 'rgba', 'ubyte', 'nearest');
+    }
+    if (!ctx.namedTextures[AlphaOrbitalsTex1]) {
+        ctx.namedTextures[AlphaOrbitalsTex1] = ctx.resources.texture('image-uint8', 'rgba', 'ubyte', 'nearest');
+    }
+
+    const values: AlphaOrbitalsSchema = {
         ...QuadValues,
-        uDimensions: ValueCell.create(params.grid.dimensions),
-        uMin: ValueCell.create(params.grid.box.min),
-        uDelta: ValueCell.create(params.grid.delta),
+        uDimensions: ValueCell.create(grid.dimensions),
+        uMin: ValueCell.create(grid.box.min),
+        uDelta: ValueCell.create(grid.delta),
         uWidth: ValueCell.create(width),
         uNCenters: ValueCell.create(data.nCenters),
         uNAlpha: ValueCell.create(data.nAlpha),
@@ -134,6 +165,9 @@ function createAlphaOrbitalsRenderable(ctx: WebGLContext, params: CollocationPar
         tCoeff: ValueCell.create({ width: data.nCoeff, height: 1, array: data.coeff }),
         tAlpha: ValueCell.create({ width: data.nAlpha, height: 1, array: data.alpha }),
         uLittleEndian: ValueCell.create(isLittleEndian()),
+        uDensity: ValueCell.create(false),
+        uOccupancy: ValueCell.create(0),
+        tCumulativeSum: ValueCell.create(ctx.namedTextures[AlphaOrbitalsTex1])
     };
 
     const schema = { ...AlphaOrbitalsSchema };
@@ -148,18 +182,18 @@ function createAlphaOrbitalsRenderable(ctx: WebGLContext, params: CollocationPar
     return createComputeRenderable(renderItem, values);
 }
 
-function getAlphaOrbitalsRenderable(ctx: WebGLContext, params: CollocationParams): AlphaOrbitalsRenderable {
+function getAlphaOrbitalsRenderable(ctx: WebGLContext, grid: CubeGridInfo, orbital: AlphaOrbital): AlphaOrbitalsRenderable {
     if (ctx.namedComputeRenderables[AlphaOrbitalsName]) {
-        const v = ctx.namedComputeRenderables[AlphaOrbitalsName].values;
+        const v = ctx.namedComputeRenderables[AlphaOrbitalsName].values as AlphaOrbitalsSchema;
 
-        const data = createTextureData(params);
+        const data = createTextureData(grid, orbital);
 
-        const [nx, ny, nz] = params.grid.dimensions;
+        const [nx, ny, nz] = grid.dimensions;
         const width = Math.ceil(Math.sqrt(nx * ny * nz));
 
-        ValueCell.update(v.uDimensions, params.grid.dimensions);
-        ValueCell.update(v.uMin, params.grid.box.min);
-        ValueCell.update(v.uDelta, params.grid.delta);
+        ValueCell.update(v.uDimensions, grid.dimensions);
+        ValueCell.update(v.uMin, grid.box.min);
+        ValueCell.update(v.uDelta, grid.delta);
         ValueCell.updateIfChanged(v.uWidth, width);
         ValueCell.updateIfChanged(v.uNCenters, data.nCenters);
         ValueCell.updateIfChanged(v.uNAlpha, data.nAlpha);
@@ -170,29 +204,25 @@ function getAlphaOrbitalsRenderable(ctx: WebGLContext, params: CollocationParams
         ValueCell.update(v.tCoeff, { width: data.nCoeff, height: 1, array: data.coeff });
         ValueCell.update(v.tAlpha, { width: data.nAlpha, height: 1, array: data.alpha });
         ValueCell.updateIfChanged(v.uLittleEndian, isLittleEndian());
+        ValueCell.updateIfChanged(v.uDensity, false);
+        ValueCell.updateIfChanged(v.uOccupancy, 0);
+        ValueCell.updateIfChanged(v.tCumulativeSum, ctx.namedTextures[AlphaOrbitalsTex1]);
 
         ctx.namedComputeRenderables[AlphaOrbitalsName].update();
     } else {
-        ctx.namedComputeRenderables[AlphaOrbitalsName] = createAlphaOrbitalsRenderable(ctx, params);
+        ctx.namedComputeRenderables[AlphaOrbitalsName] = createAlphaOrbitalsRenderable(ctx, grid, orbital);
     }
     return ctx.namedComputeRenderables[AlphaOrbitalsName];
 }
 
-export function gpuComputeAlphaOrbitalsGridValues(webgl: WebGLContext, params: CollocationParams) {
-    const [nx, ny, nz] = params.grid.dimensions;
-    const renderable = getAlphaOrbitalsRenderable(webgl, params);
+export function gpuComputeAlphaOrbitalsGridValues(webgl: WebGLContext, grid: CubeGridInfo, orbital: AlphaOrbital) {
+    const [nx, ny, nz] = grid.dimensions;
+    const renderable = getAlphaOrbitalsRenderable(webgl, grid, orbital);
     const width = renderable.values.uWidth.ref.value;
 
-    if (!webgl.namedFramebuffers[AlphaOrbitalsName]) {
-        webgl.namedFramebuffers[AlphaOrbitalsName] = webgl.resources.framebuffer();
-    }
     const framebuffer = webgl.namedFramebuffers[AlphaOrbitalsName];
-
-    if (!webgl.namedTextures[AlphaOrbitalsName]) {
-        webgl.namedTextures[AlphaOrbitalsName] = webgl.resources.texture('image-uint8', 'rgba', 'ubyte', 'nearest');
-    }
-    webgl.namedTextures[AlphaOrbitalsName].define(width, width);
-    webgl.namedTextures[AlphaOrbitalsName].attachFramebuffer(framebuffer, 'color0');
+    webgl.namedTextures[AlphaOrbitalsTex0].define(width, width);
+    webgl.namedTextures[AlphaOrbitalsTex0].attachFramebuffer(framebuffer, 'color0');
 
     const { gl, state } = webgl;
     framebuffer.bind();
@@ -211,4 +241,69 @@ export function gpuComputeAlphaOrbitalsGridValues(webgl: WebGLContext, params: C
 
 export function canComputeAlphaOrbitalsOnGPU(webgl?: WebGLContext) {
     return !!webgl?.extensions.textureFloat;
+}
+
+export async function gpuComputeAlphaOrbitalsDensityGridValues(webgl: WebGLContext, grid: CubeGridInfo, orbitals: AlphaOrbital[], ctx: RuntimeContext) {
+    await ctx.update({ message: 'Initializing...', isIndeterminate: true });
+
+    const [nx, ny, nz] = grid.dimensions;
+    const renderable = getAlphaOrbitalsRenderable(webgl, grid, orbitals[0]);
+    const width = renderable.values.uWidth.ref.value;
+
+    if (!webgl.namedFramebuffers[AlphaOrbitalsName]) {
+        webgl.namedFramebuffers[AlphaOrbitalsName] = webgl.resources.framebuffer();
+    }
+    const framebuffer = webgl.namedFramebuffers[AlphaOrbitalsName];
+    const tex = [webgl.namedTextures[AlphaOrbitalsTex0], webgl.namedTextures[AlphaOrbitalsTex1]];
+
+    tex[0].define(width, width);
+    tex[1].define(width, width);
+
+    const values = renderable.values as AlphaOrbitalsSchema;
+    const { gl, state } = webgl;
+
+    gl.viewport(0, 0, width, width);
+    gl.scissor(0, 0, width, width);
+    state.disable(gl.SCISSOR_TEST);
+    state.disable(gl.BLEND);
+    state.disable(gl.DEPTH_TEST);
+    state.depthMask(false);
+
+    gl.clearColor(0, 0, 0, 0);
+
+    tex[0].attachFramebuffer(framebuffer, 'color0');
+    gl.clear(gl.COLOR_BUFFER_BIT);
+
+    tex[1].attachFramebuffer(framebuffer, 'color0');
+    gl.clear(gl.COLOR_BUFFER_BIT);
+
+    ValueCell.update(values.uDensity, true);
+
+    const nonZero = orbitals.filter(o => o.occupancy !== 0);
+    await ctx.update({ message: 'Computing...', isIndeterminate: false, current: 0, max: nonZero.length });
+    for (let i = 0; i < nonZero.length; i++) {
+        const alpha = getNormalizedAlpha(grid.params.basis, nonZero[i].alpha, grid.params.sphericalOrder);
+
+        ValueCell.update(values.uOccupancy, nonZero[i].occupancy);
+        ValueCell.update(values.tCumulativeSum, tex[(i + 1) % 2]);
+        ValueCell.update(values.tAlpha, { width: alpha.length, height: 1, array: alpha });
+        tex[i % 2].attachFramebuffer(framebuffer, 'color0');
+        gl.viewport(0, 0, width, width);
+        gl.scissor(0, 0, width, width);
+        state.disable(gl.SCISSOR_TEST);
+        state.disable(gl.BLEND);
+        state.disable(gl.DEPTH_TEST);
+        state.depthMask(false);
+        renderable.update();
+        renderable.render();
+
+        if (i !== nonZero.length - 1 && ctx.shouldUpdate) {
+            await ctx.update({ current: i + 1 });
+        }
+    }
+
+    const array = new Uint8Array(width * width * 4);
+    webgl.readPixels(0, 0, width, width, array);
+
+    return new Float32Array(array.buffer, array.byteOffset, nx * ny * nz);
 }

src/extensions/alpha-orbitals/gpu/shader.frag.ts

@@ -28,8 +28,17 @@ uniform float uWidth;
 uniform int uNCoeff;
 uniform int uNAlpha;
 
+uniform bool uDensity;
+uniform float uOccupancy;
+uniform sampler2D tCumulativeSum;
+
 uniform bool uLittleEndian;
 
+//////////////////////////////////////////////////////////
+
+// floatToRgba adapted from https://github.com/equinor/glsl-float-to-rgba
+// MIT License, Copyright (c) 2020 Equinor
+
 float shiftRight (float v, float amt) {
   v = floor(v) + 0.5;
   return floor(v / exp2(amt));
@@ -44,8 +53,7 @@ float extractBits (float num, float from, float to) {
     from = floor(from + 0.5); to = floor(to + 0.5);
     return maskLast(shiftRight(num, from), to - from);
 }
-// Adapted from https://github.com/equinor/glsl-float-to-rgba
-// MIT License, Copyright (c) 2020 Equinor
+
 vec4 floatToRgba(float texelFloat) {
     if (texelFloat == 0.0) return vec4(0, 0, 0, 0);
     float sign = texelFloat > 0.0 ? 0.0 : 1.0;
@@ -67,6 +75,94 @@ vec4 floatToRgba(float texelFloat) {
     );
 }
 
+///////////////////////////////////////////////////////
+
+// rgbaToFloat adapted from https://github.com/ihmeuw/glsl-rgba-to-float
+// BSD 3-Clause License
+//
+// Copyright (c) 2019, Institute for Health Metrics and Evaluation All rights reserved.
+// Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+//  - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+//  - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+//  - Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
+// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+// IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
+// OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+// OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+// OF THE POSSIBILITY OF SUCH DAMAGE.
+
+ivec4 floatsToBytes(vec4 inputFloats) {
+  ivec4 bytes = ivec4(inputFloats * 255.0);
+  return (
+    uLittleEndian
+    ? bytes.abgr
+    : bytes
+  );
+}
+
+// Break the four bytes down into an array of 32 bits.
+void bytesToBits(const in ivec4 bytes, out bool bits[32]) {
+  for (int channelIndex = 0; channelIndex < 4; ++channelIndex) {
+    float acc = float(bytes[channelIndex]);
+    for (int indexInByte = 7; indexInByte >= 0; --indexInByte) {
+      float powerOfTwo = exp2(float(indexInByte));
+      bool bit = acc >= powerOfTwo;
+      bits[channelIndex * 8 + (7 - indexInByte)] = bit;
+      acc = mod(acc, powerOfTwo);
+    }
+  }
+}
+
+// Compute the exponent of the 32-bit float.
+float getExponent(bool bits[32]) {
+  const int startIndex = 1;
+  const int bitStringLength = 8;
+  const int endBeforeIndex = startIndex + bitStringLength;
+  float acc = 0.0;
+  int pow2 = bitStringLength - 1;
+  for (int bitIndex = startIndex; bitIndex < endBeforeIndex; ++bitIndex) {
+    acc += float(bits[bitIndex]) * exp2(float(pow2--));
+  }
+  return acc;
+}
+
+// Compute the mantissa of the 32-bit float.
+float getMantissa(bool bits[32], bool subnormal) {
+  const int startIndex = 9;
+  const int bitStringLength = 23;
+  const int endBeforeIndex = startIndex + bitStringLength;
+  // Leading/implicit/hidden bit convention:
+  // If the number is not subnormal (with exponent 0), we add a leading 1 digit.
+  float acc = float(!subnormal) * exp2(float(bitStringLength));
+  int pow2 = bitStringLength - 1;
+  for (int bitIndex = startIndex; bitIndex < endBeforeIndex; ++bitIndex) {
+    acc += float(bits[bitIndex]) * exp2(float(pow2--));
+  }
+  return acc;
+}
+
+// Parse the float from its 32 bits.
+float bitsToFloat(bool bits[32]) {
+  float signBit = float(bits[0]) * -2.0 + 1.0;
+  float exponent = getExponent(bits);
+  bool subnormal = abs(exponent - 0.0) < 0.01;
+  float mantissa = getMantissa(bits, subnormal);
+  float exponentBias = 127.0;
+  return signBit * mantissa * exp2(exponent - exponentBias - 23.0);
+}
+
+float rgbaToFloat(vec4 texelRGBA) {
+  ivec4 rgbaBytes = floatsToBytes(texelRGBA);
+  bool bits[32];
+  bytesToBits(rgbaBytes, bits);
+  return bitsToFloat(bits);
+}
+
+///////////////////////////////////////////////////////
+
 float L1(vec3 p, float a0, float a1, float a2) {
     return a0 * p.z + a1 * p.x + a2 * p.y;
 }
@@ -224,7 +320,12 @@ void main(void) {
         v += R(R2, coeffStart, coeffEnd, fCoeff) * Y(L, X, aO, fA);
     }
 
-    // TODO: render to single component float32 texture in WebGL2
-    gl_FragColor = floatToRgba(v);
+    
+    if (uDensity) {
+        float current = rgbaToFloat(texture2D(tCumulativeSum, gl_FragCoord.xy / vec2(uWidth, uWidth)));
+        gl_FragColor = floatToRgba(current + uOccupancy * v * v);
+    } else {
+         gl_FragColor = floatToRgba(v);
+    }
 }
 `;

src/extensions/alpha-orbitals/orbitals.ts

@@ -6,88 +6,123 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
-// gaussian:
-//     R_0, R^+_1, R^-_1, ..., R^+_l, R^-_l
-// cca:
-//     R^-_(l), R^-_(l-1), ..., R_0, ..., R^+_(l-1), R^+_l
-// cca-reverse:
-//     R^+_(l), R^+_(l-1), ..., R_0, ..., R^-_(l-1), R^-_l
-export type SphericalBasisOrder = 'gaussian' | 'cca' | 'cca-reverse';
+import { sortArray } from '../../mol-data/util';
+import { WebGLContext } from '../../mol-gl/webgl/context';
+import { Task } from '../../mol-task';
+import { sphericalCollocation } from './collocation';
+import { AlphaOrbital, createGrid, CubeGrid, CubeGridComputationParams, initCubeGrid } from './data-model';
+import { canComputeAlphaOrbitalsOnGPU, gpuComputeAlphaOrbitalsGridValues } from './gpu/compute';
 
-export function normalizeBasicOrder(
-    L: number,
-    alpha: number[],
-    order: SphericalBasisOrder
-) {
-    if (order === 'gaussian' || L === 0) return alpha;
+export function createSphericalCollocationGrid(
+    params: CubeGridComputationParams, orbital: AlphaOrbital, webgl?: WebGLContext
+): Task<CubeGrid> {
+    return Task.create('Spherical Collocation Grid', async (ctx) => {
+        const cubeGrid = initCubeGrid(params);
 
-    const ret: number[] = [alpha[L]];
-    for (let l = 0; l < L; l++) {
-        const a = alpha[L - l - 1],
-            b = alpha[L + l + 1];
-        if (order === 'cca') ret.push(b, a);
-        else ret.push(a, b);
+        let matrix: Float32Array;
+        if (canComputeAlphaOrbitalsOnGPU(webgl)) {
+            // console.time('gpu');
+            matrix = gpuComputeAlphaOrbitalsGridValues(webgl!, cubeGrid, orbital);
+            // console.timeEnd('gpu');
+        } else {
+            // console.time('cpu');
+            matrix = await sphericalCollocation(cubeGrid, orbital, ctx);
+            // console.timeEnd('cpu');
+        }
+
+        const grid = createGrid(cubeGrid, matrix, [0, 1, 2]);
+        let isovalues: { negative?: number, positive?: number } | undefined;
+
+        if (!params.doNotComputeIsovalues) {
+            isovalues = computeOrbitalIsocontourValues(matrix, 0.85);
+        }
+
+        return { grid, isovalues };
+    });
+}
+
+export function computeOrbitalIsocontourValues(input: Float32Array, cumulativeThreshold: number) {
+    let weightSum = 0;
+    for (let i = 0, _i = input.length; i < _i; i++) {
+        const v = input[i];
+        const w = v * v;
+        weightSum += w;
     }
-    return ret;
-}
+    const avgWeight = weightSum / input.length;
+    let minWeight = 3 * avgWeight;
 
-export type SphericalFunc = (
-    alpha: number[],
-    x: number,
-    y: number,
-    z: number
-) => number;
+    // do not try to identify isovalues for degenerate data
+    // e.g. all values are almost same
+    if (Math.abs(avgWeight - input[0] * input[0]) < 1e-5) {
+        return { negative: void 0, positive: void 0 };
+    }
 
-export const SphericalFunctions: SphericalFunc[] = [L0, L1, L2, L3, L4];
+    let size = 0;
+    while (true) {
+        let csum = 0;
+        size = 0;
+        for (let i = 0, _i = input.length; i < _i; i++) {
+            const v = input[i];
+            const w = v * v;
+            if (w >= minWeight) {
+                csum += w;
+                size++;
+            }
+        }
 
-// L_i functions were auto-generated.
+        if (csum / weightSum > cumulativeThreshold) {
+            break;
+        }
 
-function L0(alpha: number[], x: number, y: number, z: number) {
-    return alpha[0];
-}
+        minWeight -= avgWeight;
+    }
 
-function L1(alpha: number[], x: number, y: number, z: number) {
-    return alpha[0] * z + alpha[1] * x + alpha[2] * y;
-}
+    const values = new Float32Array(size);
+    const weights = new Float32Array(size);
+    const indices = new Int32Array(size);
 
-function L2(alpha: number[], x: number, y: number, z: number) {
-    const xx = x * x, yy = y * y, zz = z * z;
-    return (
-        alpha[0] * (-0.5 * xx - 0.5 * yy + zz) +
-        alpha[1] * (1.7320508075688772 * x * z) +
-        alpha[2] * (1.7320508075688772 * y * z) +
-        alpha[3] * (0.8660254037844386 * xx - 0.8660254037844386 * yy) +
-        alpha[4] * (1.7320508075688772 * x * y)
+    let o = 0;
+    for (let i = 0, _i = input.length; i < _i; i++) {
+        const v = input[i];
+        const w = v * v;
+        if (w >= minWeight) {
+            values[o] = v;
+            weights[o] = w;
+            indices[o] = o;
+            o++;
+        }
+    }
+
+    sortArray(
+        indices,
+        (indices, i, j) => weights[indices[j]] - weights[indices[i]]
     );
-}
 
-function L3(alpha: number[], x: number, y: number, z: number) {
-    const xx = x * x, yy = y * y, zz = z * z;
-    const xxx = xx * x, yyy = yy * y, zzz = zz * z;
-    return (
-        alpha[0] * (-1.5 * xx * z - 1.5 * yy * z + zzz) +
-        alpha[1] * (-0.6123724356957945 * xxx - 0.6123724356957945 * x * yy + 2.449489742783178 * x * zz) +
-        alpha[2] * (-0.6123724356957945 * xx * y - 0.6123724356957945 * yyy + 2.449489742783178 * y * zz) +
-        alpha[3] * (1.9364916731037085 * xx * z - 1.9364916731037085 * yy * z) +
-        alpha[4] * (3.872983346207417 * x * y * z) +
-        alpha[5] * (0.7905694150420949 * xxx - 2.3717082451262845 * x * yy) +
-        alpha[6] * (2.3717082451262845 * xx * y - 0.7905694150420949 * yyy)
-    );
-}
+    let cweight = 0,
+        cutoffIndex = 0;
+    for (let i = 0; i < size; i++) {
+        cweight += weights[indices[i]];
 
-function L4(alpha: number[], x: number, y: number, z: number) {
-    const xx = x * x, yy = y * y, zz = z * z;
-    const xxx = xx * x, yyy = yy * y, zzz = zz * z;
-    const xxxx = xxx * x, yyyy = yyy * y, zzzz = zzz * z;
-    return (
-        alpha[0] * (0.375 * xxxx + 0.75 * xx * yy + 0.375 * yyyy - 3.0 * xx * zz - 3.0 * yy * zz + zzzz) +
-        alpha[1] * (-2.3717082451262845 * xxx * z - 2.3717082451262845 * x * yy * z + 3.1622776601683795 * x * zzz) +
-        alpha[2] * (-2.3717082451262845 * xx * y * z - 2.3717082451262845 * yyy * z + 3.1622776601683795 * y * zzz) +
-        alpha[3] * (-0.5590169943749475 * xxxx + 0.5590169943749475 * yyyy + 3.3541019662496847 * xx * zz - 3.3541019662496847 * yy * zz) +
-        alpha[4] * (-1.118033988749895 * xxx * y - 1.118033988749895 * x * yyy + 6.708203932499369 * x * y * zz) +
-        alpha[5] * (2.091650066335189 * xxx * z + -6.274950199005566 * x * yy * z) +
-        alpha[6] * (6.274950199005566 * xx * y * z + -2.091650066335189 * yyy * z) +
-        alpha[7] * (0.739509972887452 * xxxx - 4.437059837324712 * xx * yy + 0.739509972887452 * yyyy) +
-        alpha[8] * (2.958039891549808 * xxx * y + -2.958039891549808 * x * yyy)
-    );
-}
+        if (cweight / weightSum >= cumulativeThreshold) {
+            cutoffIndex = i;
+            break;
+        }
+    }
+
+    let positive = Number.POSITIVE_INFINITY,
+        negative = Number.NEGATIVE_INFINITY;
+
+    for (let i = 0; i < cutoffIndex; i++) {
+        const v = values[indices[i]];
+        if (v > 0) {
+            if (v < positive) positive = v;
+        } else if (v < 0) {
+            if (v > negative) negative = v;
+        }
+    }
+
+    return {
+        negative: negative !== Number.NEGATIVE_INFINITY ? negative : void 0,
+        positive: positive !== Number.POSITIVE_INFINITY ? positive : void 0,
+    };
+}

src/extensions/alpha-orbitals/spherical-functions.ts

@@ -0,0 +1,93 @@
+/**
+ * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * Inspired by https://github.com/dgasmith/gau2grid.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+// gaussian:
+//     R_0, R^+_1, R^-_1, ..., R^+_l, R^-_l
+// cca:
+//     R^-_(l), R^-_(l-1), ..., R_0, ..., R^+_(l-1), R^+_l
+// cca-reverse:
+//     R^+_(l), R^+_(l-1), ..., R_0, ..., R^-_(l-1), R^-_l
+export type SphericalBasisOrder = 'gaussian' | 'cca' | 'cca-reverse';
+
+export function normalizeBasicOrder(
+    L: number,
+    alpha: number[],
+    order: SphericalBasisOrder
+) {
+    if (order === 'gaussian' || L === 0) return alpha;
+
+    const ret: number[] = [alpha[L]];
+    for (let l = 0; l < L; l++) {
+        const a = alpha[L - l - 1],
+            b = alpha[L + l + 1];
+        if (order === 'cca') ret.push(b, a);
+        else ret.push(a, b);
+    }
+    return ret;
+}
+
+export type SphericalFunc = (
+    alpha: number[],
+    x: number,
+    y: number,
+    z: number
+) => number;
+
+export const SphericalFunctions: SphericalFunc[] = [L0, L1, L2, L3, L4];
+
+// L_i functions were auto-generated.
+
+function L0(alpha: number[], x: number, y: number, z: number) {
+    return alpha[0];
+}
+
+function L1(alpha: number[], x: number, y: number, z: number) {
+    return alpha[0] * z + alpha[1] * x + alpha[2] * y;
+}
+
+function L2(alpha: number[], x: number, y: number, z: number) {
+    const xx = x * x, yy = y * y, zz = z * z;
+    return (
+        alpha[0] * (-0.5 * xx - 0.5 * yy + zz) +
+        alpha[1] * (1.7320508075688772 * x * z) +
+        alpha[2] * (1.7320508075688772 * y * z) +
+        alpha[3] * (0.8660254037844386 * xx - 0.8660254037844386 * yy) +
+        alpha[4] * (1.7320508075688772 * x * y)
+    );
+}
+
+function L3(alpha: number[], x: number, y: number, z: number) {
+    const xx = x * x, yy = y * y, zz = z * z;
+    const xxx = xx * x, yyy = yy * y, zzz = zz * z;
+    return (
+        alpha[0] * (-1.5 * xx * z - 1.5 * yy * z + zzz) +
+        alpha[1] * (-0.6123724356957945 * xxx - 0.6123724356957945 * x * yy + 2.449489742783178 * x * zz) +
+        alpha[2] * (-0.6123724356957945 * xx * y - 0.6123724356957945 * yyy + 2.449489742783178 * y * zz) +
+        alpha[3] * (1.9364916731037085 * xx * z - 1.9364916731037085 * yy * z) +
+        alpha[4] * (3.872983346207417 * x * y * z) +
+        alpha[5] * (0.7905694150420949 * xxx - 2.3717082451262845 * x * yy) +
+        alpha[6] * (2.3717082451262845 * xx * y - 0.7905694150420949 * yyy)
+    );
+}
+
+function L4(alpha: number[], x: number, y: number, z: number) {
+    const xx = x * x, yy = y * y, zz = z * z;
+    const xxx = xx * x, yyy = yy * y, zzz = zz * z;
+    const xxxx = xxx * x, yyyy = yyy * y, zzzz = zzz * z;
+    return (
+        alpha[0] * (0.375 * xxxx + 0.75 * xx * yy + 0.375 * yyyy - 3.0 * xx * zz - 3.0 * yy * zz + zzzz) +
+        alpha[1] * (-2.3717082451262845 * xxx * z - 2.3717082451262845 * x * yy * z + 3.1622776601683795 * x * zzz) +
+        alpha[2] * (-2.3717082451262845 * xx * y * z - 2.3717082451262845 * yyy * z + 3.1622776601683795 * y * zzz) +
+        alpha[3] * (-0.5590169943749475 * xxxx + 0.5590169943749475 * yyyy + 3.3541019662496847 * xx * zz - 3.3541019662496847 * yy * zz) +
+        alpha[4] * (-1.118033988749895 * xxx * y - 1.118033988749895 * x * yyy + 6.708203932499369 * x * y * zz) +
+        alpha[5] * (2.091650066335189 * xxx * z + -6.274950199005566 * x * yy * z) +
+        alpha[6] * (6.274950199005566 * xx * y * z + -2.091650066335189 * yyy * z) +
+        alpha[7] * (0.739509972887452 * xxxx - 4.437059837324712 * xx * yy + 0.739509972887452 * yyyy) +
+        alpha[8] * (2.958039891549808 * xxx * y + -2.958039891549808 * x * yyy)
+    );
+}

src/extensions/alpha-orbitals/transforms.ts

@@ -5,11 +5,11 @@
  */
 
 import { PluginStateObject, PluginStateTransform } from '../../mol-plugin-state/objects';
-import { Basis, createSphericalCollocationGrid, CubeGrid } from './cubes';
+import { createSphericalCollocationGrid } from './orbitals';
 import { ParamDefinition as PD } from '../../mol-util/param-definition';
 import { Task } from '../../mol-task';
 import { CustomProperties } from '../../mol-model/custom-property';
-import { SphericalBasisOrder } from './orbitals';
+import { SphericalBasisOrder } from './spherical-functions';
 import { Volume } from '../../mol-model/volume';
 import { PluginContext } from '../../mol-plugin/context';
 import { ColorNames } from '../../mol-util/color/names';
@@ -17,8 +17,10 @@ import { createVolumeRepresentationParams } from '../../mol-plugin-state/helpers
 import { StateTransformer } from '../../mol-state';
 import { Theme } from '../../mol-theme/theme';
 import { VolumeRepresentation3DHelpers } from '../../mol-plugin-state/transforms/representation';
+import { AlphaOrbital, Basis, CubeGrid } from './data-model';
+import { createSphericalCollocationDensityGrid } from './density';
 
-export class BasisAndOrbitals extends PluginStateObject.Create<{ basis: Basis, order: SphericalBasisOrder, orbitals: { energy: number, alpha: number[] }[] }>({ name: 'Basis', typeClass: 'Object' }) { }
+export class BasisAndOrbitals extends PluginStateObject.Create<{ basis: Basis, order: SphericalBasisOrder, orbitals: AlphaOrbital[] }>({ name: 'Basis', typeClass: 'Object' }) { }
 
 export const StaticBasisAndOrbitals = PluginStateTransform.BuiltIn({
     name: 'static-basis-and-orbitals',
@@ -29,7 +31,7 @@ export const StaticBasisAndOrbitals = PluginStateTransform.BuiltIn({
         label: PD.Text('Orbital Data', { isHidden: true }),
         basis: PD.Value<Basis>(void 0 as any, { isHidden: true }),
         order: PD.Text<SphericalBasisOrder>('gaussian' as SphericalBasisOrder, { isHidden: true }),
-        orbitals: PD.Value<{ energy: number, alpha: number[] }[]>([], { isHidden: true })
+        orbitals: PD.Value<AlphaOrbital[]>([], { isHidden: true })
     },
 })({
     apply({ params }) {
@@ -38,7 +40,6 @@ export const StaticBasisAndOrbitals = PluginStateTransform.BuiltIn({
 });
 
 const CreateOrbitalVolumeParamBase = {
-    forceCpuCompute: PD.Boolean(false),
     cutoffThreshold: PD.Numeric(0.0015, { min: 0, max: 0.1, step: 0.0001 }),
     boxExpand: PD.Numeric(4.5, { min: 0, max: 7, step: 0.1 }),
     gridSpacing: PD.ObjectList({ atomCount: PD.Numeric(0), spacing: PD.Numeric(0.35, { min: 0.1, max: 2, step: 0.01 }) }, e => `Atoms ${e.atomCount}: ${e.spacing}`, {
@@ -72,11 +73,38 @@ export const CreateOrbitalVolume = PluginStateTransform.BuiltIn({
             const data = await createSphericalCollocationGrid({
                 basis: a.data.basis,
                 cutoffThreshold: params.cutoffThreshold,
-                alphaOrbitals: a.data.orbitals[params.index].alpha,
                 sphericalOrder: a.data.order,
                 boxExpand: params.boxExpand,
                 gridSpacing: params.gridSpacing.map(e => [e.atomCount, e.spacing] as [number, number])
-            }, params.forceCpuCompute ? void 0 : plugin.canvas3d?.webgl).runInContext(ctx);
+            }, a.data.orbitals[params.index], plugin.canvas3d?.webgl).runInContext(ctx);
+            const volume: Volume = {
+                grid: data.grid,
+                sourceData: { name: 'custom grid', kind: 'alpha-orbitals', data },
+                customProperties: new CustomProperties(),
+                _propertyData: Object.create(null),
+            };
+
+            return new PluginStateObject.Volume.Data(volume, { label: 'Orbital Volume' });
+        });
+    }
+});
+
+export const CreateOrbitalDensityVolume = PluginStateTransform.BuiltIn({
+    name: 'create-orbital-density-volume',
+    display: 'Orbital Density Volume',
+    from: BasisAndOrbitals,
+    to: PluginStateObject.Volume.Data,
+    params: CreateOrbitalVolumeParamBase
+})({
+    apply({ a, params }, plugin: PluginContext) {
+        return Task.create('Orbital Volume', async ctx => {
+            const data = await createSphericalCollocationDensityGrid({
+                basis: a.data.basis,
+                cutoffThreshold: params.cutoffThreshold,
+                sphericalOrder: a.data.order,
+                boxExpand: params.boxExpand,
+                gridSpacing: params.gridSpacing.map(e => [e.atomCount, e.spacing] as [number, number])
+            }, a.data.orbitals, plugin.canvas3d?.webgl).runInContext(ctx);
             const volume: Volume = {
                 grid: data.grid,
                 sourceData: { name: 'custom grid', kind: 'alpha-orbitals', data },

src/extensions/rcsb/validation-report/color/geometry-quality.ts

@@ -32,7 +32,7 @@ const ColorLegend = TableLegend([
 ]);
 
 export function getGeometricQualityColorThemeParams(ctx: ThemeDataContext) {
-    const validationReport = ctx.structure && ValidationReportProvider.get(ctx.structure.models[0]).value;
+    const validationReport = !!ctx.structure && ctx.structure.models.length > 0 && ValidationReportProvider.get(ctx.structure.models[0]).value;
     const options: [string, string][] = [];
     if (validationReport) {
         const kinds = new Set<string>();
@@ -48,7 +48,7 @@ export type GeometricQualityColorThemeParams = ReturnType<typeof getGeometricQua
 export function GeometryQualityColorTheme(ctx: ThemeDataContext, props: PD.Values<GeometricQualityColorThemeParams>): ColorTheme<GeometricQualityColorThemeParams> {
     let color: LocationColor = () => DefaultColor;
 
-    const validationReport = ctx.structure && ValidationReportProvider.get(ctx.structure.models[0]);
+    const validationReport = !!ctx.structure && ctx.structure.models.length > 0 ? ValidationReportProvider.get(ctx.structure.models[0]) : void 0;
     const contextHash = validationReport?.version;
 
     const value = validationReport?.value;

src/mol-canvas3d/canvas3d.ts

@@ -149,10 +149,10 @@ namespace Canvas3D {
     export interface DragEvent { current: Representation.Loci, buttons: ButtonsType, button: ButtonsType.Flag, modifiers: ModifiersKeys, pageStart: Vec2, pageEnd: Vec2 }
     export interface ClickEvent { current: Representation.Loci, buttons: ButtonsType, button: ButtonsType.Flag, modifiers: ModifiersKeys, page?: Vec2, position?: Vec3 }
 
-    export function fromCanvas(canvas: HTMLCanvasElement, props: Partial<Canvas3DProps> = {}, attribs: Partial<{ antialias: boolean, pixelScale?: number, pickScale: number, enableWboit: boolean }> = {}) {
+    export function fromCanvas(canvas: HTMLCanvasElement, props: Partial<Canvas3DProps> = {}, attribs: Partial<{ antialias: boolean, pixelScale: number, pickScale: number, enableWboit: boolean }> = {}) {
         const gl = getGLContext(canvas, {
             alpha: true,
-            antialias: attribs.antialias ?? true,
+            antialias: (attribs.antialias ?? true) && !attribs.enableWboit,
             depth: true,
             preserveDrawingBuffer: true,
             premultipliedAlpha: true,
@@ -420,7 +420,8 @@ namespace Canvas3D {
                 const b = r.values.boundingSphere.ref.value;
                 if (!b.radius) continue;
 
-                if (!Sphere3D.includes(oldBoundingSphereVisible, b)) return true;
+                const cameraDist = Vec3.distance(cameraSphere.center, b.center);
+                if ((cameraDist > cameraSphere.radius || cameraDist > b.radius || b.radius > camera.state.radiusMax) && !Sphere3D.includes(oldBoundingSphereVisible, b)) return true;
                 if (Sphere3D.overlaps(cameraSphere, b)) cameraSphereOverlapsNone = false;
             }
 

src/mol-canvas3d/helper/interaction-events.ts

@@ -93,7 +93,7 @@ export class Canvas3dInteractionHelper {
 
     private leave() {
         this.inside = false;
-        if (Representation.Loci.isEmpty(this.prevLoci)) {
+        if (!Representation.Loci.isEmpty(this.prevLoci)) {
             this.prevLoci = Representation.Loci.Empty;
             this.events.hover.next({ current: this.prevLoci, buttons: this.buttons, button: this.button, modifiers: this.modifiers });
         }

src/mol-canvas3d/passes/draw.ts

@@ -75,7 +75,7 @@ export class DrawPass {
 
         this.drawTarget = createNullRenderTarget(webgl.gl);
 
-        this.colorTarget = webgl.createRenderTarget(width, height);
+        this.colorTarget = webgl.createRenderTarget(width, height, true, 'uint8', 'linear');
         this.packedDepth = !extensions.depthTexture;
 
         this.depthTarget = webgl.createRenderTarget(width, height);
@@ -195,7 +195,7 @@ export class DrawPass {
         if (!toDrawingBuffer) {
             if (!this.packedDepth) {
                 this.depthTextureVolumes.attachFramebuffer(this.colorTarget.framebuffer, 'depth');
-                renderer.clearDepth();
+                renderer.clearDepth(); // from previous frame
             }
             renderer.renderBlendedVolume(scene.volumes, camera, this.depthTexturePrimitives);
 
@@ -206,6 +206,10 @@ export class DrawPass {
                 renderer.renderDepth(scene.volumes, camera, this.depthTexturePrimitives);
                 this.colorTarget.bind();
             }
+
+            if (!this.packedDepth) {
+                this.depthTexturePrimitives.attachFramebuffer(this.colorTarget.framebuffer, 'depth');
+            }
         }
 
         renderer.renderBlendedTransparent(scene.primitives, camera, null);

src/mol-canvas3d/passes/multi-sample.ts

@@ -193,7 +193,7 @@ export class MultiSamplePass {
 
         const { x, y, width, height } = camera.viewport;
         const sampleWeight = 1.0 / offsetList.length;
-        const postprocessingEnabled = PostprocessingPass.isEnabled(props.postprocessing);
+        const postprocessingEnabled = PostprocessingPass.isEnabled(props.postprocessing) || props.postprocessing.antialiasing.name === 'on';
 
         if (sampleIndex === -1) {
             drawPass.render(renderer, camera, scene, helper, false, transparentBackground);

src/mol-canvas3d/passes/postprocessing.ts

@@ -17,9 +17,9 @@ import { ParamDefinition as PD } from '../../mol-util/param-definition';
 import { RenderTarget } from '../../mol-gl/webgl/render-target';
 import { DrawPass } from './draw';
 import { Camera, ICamera } from '../../mol-canvas3d/camera';
-
 import quad_vert from '../../mol-gl/shader/quad.vert';
 import postprocessing_frag from '../../mol-gl/shader/postprocessing.frag';
+import fxaa_frag from '../../mol-gl/shader/fxaa.frag';
 import { StereoCamera } from '../camera/stereo';
 
 const PostprocessingSchema = {
@@ -48,11 +48,14 @@ const PostprocessingShaderCode = ShaderCode('postprocessing', quad_vert, postpro
 type PostprocessingRenderable = ComputeRenderable<Values<typeof PostprocessingSchema>>
 
 function getPostprocessingRenderable(ctx: WebGLContext, colorTexture: Texture, depthTexture: Texture): PostprocessingRenderable {
+    const width = colorTexture.getWidth();
+    const height = colorTexture.getHeight();
+
     const values: Values<typeof PostprocessingSchema> = {
         ...QuadValues,
         tColor: ValueCell.create(colorTexture),
         tPackedDepth: ValueCell.create(depthTexture),
-        uTexSize: ValueCell.create(Vec2.create(colorTexture.getWidth(), colorTexture.getHeight())),
+        uTexSize: ValueCell.create(Vec2.create(width, height)),
 
         dOrthographic: ValueCell.create(0),
         uNear: ValueCell.create(1),
@@ -92,7 +95,16 @@ export const PostprocessingParams = {
             threshold: PD.Numeric(0.8, { min: 0, max: 5, step: 0.01 }),
         }),
         off: PD.Group({})
-    }, { cycle: true, description: 'Draw outline around 3D objects' })
+    }, { cycle: true, description: 'Draw outline around 3D objects' }),
+    antialiasing: PD.MappedStatic('on', {
+        on: PD.Group({
+            edgeThresholdMin:PD.Numeric(0.0312, { min: 0.0312, max: 0.0833, step: 0.0001 }, { description: 'Trims the algorithm from processing darks.' }),
+            edgeThresholdMax: PD.Numeric(0.063, { min: 0.063, max: 0.333, step: 0.001 }, { description: 'The minimum amount of local contrast required to apply algorithm.' }),
+            iterations: PD.Numeric(12, { min: 0, max: 32, step: 1 }, { description: 'Number of edge exploration steps.' }),
+            subpixelQuality: PD.Numeric(1.00, { min: 0.00, max: 1.00, step: 0.01 }, { description: 'Choose the amount of sub-pixel aliasing removal.' }),
+        }),
+        off: PD.Group({})
+    }, { cycle: true, description: 'Fast Approximate Anti-Aliasing (FXAA)' }),
 };
 export type PostprocessingProps = PD.Values<typeof PostprocessingParams>
 
@@ -101,21 +113,21 @@ export class PostprocessingPass {
         return props.occlusion.name === 'on' || props.outline.name === 'on';
     }
 
-    target: RenderTarget
-    renderable: PostprocessingRenderable
+    readonly target: RenderTarget
+
+    private readonly tmpTarget: RenderTarget
+    private readonly renderable: PostprocessingRenderable
+    private readonly fxaa: FxaaRenderable
 
     constructor(private webgl: WebGLContext, private drawPass: DrawPass) {
-        this.target = webgl.createRenderTarget(drawPass.colorTarget.getWidth(), drawPass.colorTarget.getHeight(), false);
         const { colorTarget, depthTexture } = drawPass;
-        this.renderable = getPostprocessingRenderable(webgl, colorTarget.texture, depthTexture);
-    }
+        const width = colorTarget.getWidth();
+        const height = colorTarget.getHeight();
 
-    private bindTarget(toDrawingBuffer: boolean) {
-        if (toDrawingBuffer) {
-            this.webgl.unbindFramebuffer();
-        } else {
-            this.target.bind();
-        }
+        this.target = webgl.createRenderTarget(width, height, false);
+        this.tmpTarget = webgl.createRenderTarget(width, height, false, 'uint8', 'linear');
+        this.renderable = getPostprocessingRenderable(webgl, colorTarget.texture, depthTexture);
+        this.fxaa = getFxaaRenderable(webgl, this.tmpTarget.texture);
     }
 
     syncSize() {
@@ -125,11 +137,29 @@ export class PostprocessingPass {
         const [w, h] = this.renderable.values.uTexSize.ref.value;
         if (width !== w || height !== h) {
             this.target.setSize(width, height);
+            this.tmpTarget.setSize(width, height);
             ValueCell.update(this.renderable.values.uTexSize, Vec2.set(this.renderable.values.uTexSize.ref.value, width, height));
+            ValueCell.update(this.fxaa.values.uTexSizeInv, Vec2.set(this.fxaa.values.uTexSizeInv.ref.value, 1 / width, 1 / height));
         }
     }
 
-    _render(camera: ICamera, toDrawingBuffer: boolean, props: PostprocessingProps) {
+    private updateState(camera: ICamera) {
+        const { gl, state } = this.webgl;
+
+        state.disable(gl.SCISSOR_TEST);
+        state.disable(gl.BLEND);
+        state.disable(gl.DEPTH_TEST);
+        state.depthMask(false);
+
+        const { x, y, width, height } = camera.viewport;
+        gl.viewport(x, y, width, height);
+        gl.scissor(x, y, width, height);
+
+        state.clearColor(0, 0, 0, 1);
+        gl.clear(gl.COLOR_BUFFER_BIT);
+    }
+
+    private _renderPostprocessing(camera: ICamera, toDrawingBuffer: boolean, props: PostprocessingProps) {
         const { values } = this.renderable;
 
         ValueCell.updateIfChanged(values.uFar, camera.far);
@@ -166,21 +196,66 @@ export class PostprocessingPass {
             this.renderable.update();
         }
 
-        const { gl, state } = this.webgl;
-        this.bindTarget(toDrawingBuffer);
-
-        state.disable(gl.SCISSOR_TEST);
-        state.disable(gl.BLEND);
-        state.disable(gl.DEPTH_TEST);
-        state.depthMask(false);
-
-        const { x, y, width, height } = camera.viewport;
-        gl.viewport(x, y, width, height);
-        gl.scissor(x, y, width, height);
+        if (props.antialiasing.name === 'on') {
+            this.tmpTarget.bind();
+        } else if (toDrawingBuffer) {
+            this.webgl.unbindFramebuffer();
+        } else {
+            this.target.bind();
+        }
 
+        this.updateState(camera);
         this.renderable.render();
     }
 
+    private _renderFxaa(camera: ICamera, toDrawingBuffer: boolean, props: PostprocessingProps) {
+        if (props.antialiasing.name === 'off') return;
+
+        const { values } = this.fxaa;
+
+        let needsUpdate = false;
+
+        const input = (props.occlusion.name === 'on' || props.outline.name === 'on')
+            ? this.tmpTarget.texture : this.drawPass.colorTarget.texture;
+        if (values.tColor.ref.value !== input) {
+            ValueCell.update(this.fxaa.values.tColor, input);
+            needsUpdate = true;
+        }
+
+        const { edgeThresholdMin, edgeThresholdMax, iterations, subpixelQuality } = props.antialiasing.params;
+        if (values.dEdgeThresholdMin.ref.value !== edgeThresholdMin) needsUpdate = true;
+        ValueCell.updateIfChanged(values.dEdgeThresholdMin, edgeThresholdMin);
+        if (values.dEdgeThresholdMax.ref.value !== edgeThresholdMax) needsUpdate = true;
+        ValueCell.updateIfChanged(values.dEdgeThresholdMax, edgeThresholdMax);
+        if (values.dIterations.ref.value !== iterations) needsUpdate = true;
+        ValueCell.updateIfChanged(values.dIterations, iterations);
+        if (values.dSubpixelQuality.ref.value !== subpixelQuality) needsUpdate = true;
+        ValueCell.updateIfChanged(values.dSubpixelQuality, subpixelQuality);
+
+        if (needsUpdate) {
+            this.fxaa.update();
+        }
+
+        if (toDrawingBuffer) {
+            this.webgl.unbindFramebuffer();
+        } else {
+            this.target.bind();
+        }
+
+        this.updateState(camera);
+        this.fxaa.render();
+    }
+
+    private _render(camera: ICamera, toDrawingBuffer: boolean, props: PostprocessingProps) {
+        if (props.occlusion.name === 'on' || props.outline.name === 'on' || props.antialiasing.name === 'off') {
+            this._renderPostprocessing(camera, toDrawingBuffer, props);
+        }
+
+        if (props.antialiasing.name === 'on') {
+            this._renderFxaa(camera, toDrawingBuffer, props);
+        }
+    }
+
     render(camera: Camera | StereoCamera, toDrawingBuffer: boolean, props: PostprocessingProps) {
         if (StereoCamera.is(camera)) {
             this._render(camera.left, toDrawingBuffer, props);
@@ -189,4 +264,40 @@ export class PostprocessingPass {
             this._render(camera, toDrawingBuffer, props);
         }
     }
-}
+}
+
+//
+
+const FxaaSchema = {
+    ...QuadSchema,
+    tColor: TextureSpec('texture', 'rgba', 'ubyte', 'nearest'),
+    uTexSizeInv: UniformSpec('v2'),
+
+    dEdgeThresholdMin: DefineSpec('number'),
+    dEdgeThresholdMax: DefineSpec('number'),
+    dIterations: DefineSpec('number'),
+    dSubpixelQuality: DefineSpec('number'),
+};
+const FxaaShaderCode = ShaderCode('fxaa', quad_vert, fxaa_frag);
+type FxaaRenderable = ComputeRenderable<Values<typeof FxaaSchema>>
+
+function getFxaaRenderable(ctx: WebGLContext, colorTexture: Texture): FxaaRenderable {
+    const width = colorTexture.getWidth();
+    const height = colorTexture.getHeight();
+
+    const values: Values<typeof FxaaSchema> = {
+        ...QuadValues,
+        tColor: ValueCell.create(colorTexture),
+        uTexSizeInv: ValueCell.create(Vec2.create(1 / width, 1 / height)),
+
+        dEdgeThresholdMin: ValueCell.create(0.0312),
+        dEdgeThresholdMax: ValueCell.create(0.125),
+        dIterations: ValueCell.create(12),
+        dSubpixelQuality: ValueCell.create(0.75),
+    };
+
+    const schema = { ...FxaaSchema };
+    const renderItem = createComputeRenderItem(ctx, 'triangles', FxaaShaderCode, schema, values);
+
+    return createComputeRenderable(renderItem, values);
+}

src/mol-gl/renderer.ts

@@ -422,7 +422,6 @@ namespace Renderer {
 
         const renderBlendedTransparent = (group: Scene.Group, camera: ICamera, depthTexture: Texture | null) => {
             state.enable(gl.DEPTH_TEST);
-            state.depthMask(true);
 
             updateInternal(group, camera, depthTexture, false);
 

src/mol-gl/shader-code.ts

@@ -151,7 +151,7 @@ export type ShaderDefines = {
     [k: string]: ValueCell<DefineType>
 }
 
-function getDefinesCode (defines: ShaderDefines) {
+function getDefinesCode(defines: ShaderDefines) {
     if (defines === undefined) return '';
     const lines = [];
     for (const name in defines) {

src/mol-gl/shader/chunks/wboit-write.glsl.ts

@@ -5,7 +5,8 @@ export default `
             discard;
         }
     } else if (uRenderWboit) {
-        if (preFogAlpha != 1.0 && !interior && fragmentDepth < getDepth(gl_FragCoord.xy / uDrawingBufferSize)) {
+        // the 'fragmentDepth > 0.99' check is to handle precision issues with packed depth
+        if (preFogAlpha != 1.0 && !interior && (fragmentDepth < getDepth(gl_FragCoord.xy / uDrawingBufferSize) || fragmentDepth > 0.99)) {
             float alpha = gl_FragColor.a;
             float wboitWeight = alpha * clamp(pow(1.0 - fragmentDepth, 2.0), 0.01, 1.0);
             gl_FragColor = vec4(gl_FragColor.rgb * alpha * wboitWeight, alpha);

src/mol-gl/shader/fxaa.frag.ts

@@ -0,0 +1,229 @@
+export default `
+precision highp float;
+precision highp int;
+precision highp sampler2D;
+
+uniform sampler2D tColor;
+uniform vec2 uTexSizeInv;
+
+// adapted from https://github.com/kosua20/Rendu
+// MIT License Copyright (c) 2017 Simon Rodriguez
+
+#define QUALITY(q) ((q) < 5 ? 1.0 : ((q) > 5 ? ((q) < 10 ? 2.0 : ((q) < 11 ? 4.0 : 8.0)) : 1.5))
+
+float rgb2luma(vec3 rgb){
+    return sqrt(dot(rgb, vec3(0.299, 0.587, 0.114)));
+}
+
+float sampleLuma(vec2 uv) {
+    return rgb2luma(texture2D(tColor, uv).rgb);
+}
+
+float sampleLuma(vec2 uv, float uOffset, float vOffset) {
+    uv += uTexSizeInv * vec2(uOffset, vOffset);
+    return sampleLuma(uv);
+}
+
+void main(void) {
+    vec2 coords = gl_FragCoord.xy * uTexSizeInv;
+    vec2 inverseScreenSize = uTexSizeInv;
+
+    vec4 colorCenter = texture2D(tColor, coords);
+
+    // Luma at the current fragment
+    float lumaCenter = rgb2luma(colorCenter.rgb);
+
+    // Luma at the four direct neighbours of the current fragment.
+    float lumaDown = sampleLuma(coords, 0.0, -1.0);
+    float lumaUp = sampleLuma(coords, 0.0, 1.0);
+    float lumaLeft = sampleLuma(coords, -1.0, 0.0);
+    float lumaRight = sampleLuma(coords, 1.0, 0.0);
+
+    // Find the maximum and minimum luma around the current fragment.
+    float lumaMin = min(lumaCenter, min(min(lumaDown, lumaUp), min(lumaLeft, lumaRight)));
+    float lumaMax = max(lumaCenter, max(max(lumaDown, lumaUp), max(lumaLeft, lumaRight)));
+
+    // Compute the delta.
+    float lumaRange = lumaMax - lumaMin;
+
+    // If the luma variation is lower that a threshold (or if we are in a really dark area),
+    // we are not on an edge, don't perform any AA.
+    if (lumaRange < max(dEdgeThresholdMin, lumaMax * dEdgeThresholdMax)) {
+        gl_FragColor = colorCenter;
+        return;
+    }
+
+    // Query the 4 remaining corners lumas.
+    float lumaDownLeft = sampleLuma(coords, -1.0, -1.0);
+    float lumaUpRight = sampleLuma(coords, 1.0, 1.0);
+    float lumaUpLeft = sampleLuma(coords, -1.0, 1.0);
+    float lumaDownRight = sampleLuma(coords, 1.0, -1.0);
+
+    // Combine the four edges lumas (using intermediary variables for future computations
+    // with the same values).
+    float lumaDownUp = lumaDown + lumaUp;
+    float lumaLeftRight = lumaLeft + lumaRight;
+
+    // Same for corners
+    float lumaLeftCorners = lumaDownLeft + lumaUpLeft;
+    float lumaDownCorners = lumaDownLeft + lumaDownRight;
+    float lumaRightCorners = lumaDownRight + lumaUpRight;
+    float lumaUpCorners = lumaUpRight + lumaUpLeft;
+
+    // Compute an estimation of the gradient along the horizontal and vertical axis.
+    float edgeHorizontal = abs(-2.0 * lumaLeft + lumaLeftCorners) + abs(-2.0 * lumaCenter + lumaDownUp) * 2.0 + abs(-2.0 * lumaRight + lumaRightCorners);
+    float edgeVertical = abs(-2.0 * lumaUp + lumaUpCorners) + abs(-2.0 * lumaCenter + lumaLeftRight) * 2.0 + abs(-2.0 * lumaDown + lumaDownCorners);
+
+    // Is the local edge horizontal or vertical ?
+    bool isHorizontal = (edgeHorizontal >= edgeVertical);
+
+    // Choose the step size (one pixel) accordingly.
+    float stepLength = isHorizontal ? inverseScreenSize.y : inverseScreenSize.x;
+
+    // Select the two neighboring texels lumas in the opposite direction to the local edge.
+    float luma1 = isHorizontal ? lumaDown : lumaLeft;
+    float luma2 = isHorizontal ? lumaUp : lumaRight;
+    // Compute gradients in this direction.
+    float gradient1 = luma1 - lumaCenter;
+    float gradient2 = luma2 - lumaCenter;
+
+    // Which direction is the steepest ?
+    bool is1Steepest = abs(gradient1) >= abs(gradient2);
+
+    // Gradient in the corresponding direction, normalized.
+    float gradientScaled = 0.25 * max(abs(gradient1), abs(gradient2));
+
+    // Average luma in the correct direction.
+    float lumaLocalAverage = 0.0;
+    if(is1Steepest){
+        // Switch the direction
+        stepLength = -stepLength;
+        lumaLocalAverage = 0.5 * (luma1 + lumaCenter);
+    } else {
+        lumaLocalAverage = 0.5 * (luma2 + lumaCenter);
+    }
+
+    // Shift UV in the correct direction by half a pixel.
+    vec2 currentUv = coords;
+    if(isHorizontal){
+        currentUv.y += stepLength * 0.5;
+    } else {
+        currentUv.x += stepLength * 0.5;
+    }
+
+    // Compute offset (for each iteration step) in the right direction.
+    vec2 offset = isHorizontal ? vec2(inverseScreenSize.x, 0.0) : vec2(0.0, inverseScreenSize.y);
+    // Compute UVs to explore on each side of the edge, orthogonally.
+    // The QUALITY allows us to step faster.
+    vec2 uv1 = currentUv - offset * QUALITY(0);
+    vec2 uv2 = currentUv + offset * QUALITY(0);
+
+    // Read the lumas at both current extremities of the exploration segment,
+    // and compute the delta wrt to the local average luma.
+    float lumaEnd1 = sampleLuma(uv1);
+    float lumaEnd2 = sampleLuma(uv2);
+    lumaEnd1 -= lumaLocalAverage;
+    lumaEnd2 -= lumaLocalAverage;
+
+    // If the luma deltas at the current extremities is larger than the local gradient,
+    // we have reached the side of the edge.
+    bool reached1 = abs(lumaEnd1) >= gradientScaled;
+    bool reached2 = abs(lumaEnd2) >= gradientScaled;
+    bool reachedBoth = reached1 && reached2;
+
+    // If the side is not reached, we continue to explore in this direction.
+    if(!reached1){
+        uv1 -= offset * QUALITY(1);
+    }
+    if(!reached2){
+        uv2 += offset * QUALITY(1);
+    }
+
+    // If both sides have not been reached, continue to explore.
+    if(!reachedBoth){
+        for(int i = 2; i < dIterations; i++){
+            // If needed, read luma in 1st direction, compute delta.
+            if(!reached1){
+                lumaEnd1 = sampleLuma(uv1);
+                lumaEnd1 = lumaEnd1 - lumaLocalAverage;
+            }
+            // If needed, read luma in opposite direction, compute delta.
+            if(!reached2){
+                lumaEnd2 = sampleLuma(uv2);
+                lumaEnd2 = lumaEnd2 - lumaLocalAverage;
+            }
+            // If the luma deltas at the current extremities is larger than the local gradient,
+            // we have reached the side of the edge.
+            reached1 = abs(lumaEnd1) >= gradientScaled;
+            reached2 = abs(lumaEnd2) >= gradientScaled;
+            reachedBoth = reached1 && reached2;
+
+            // If the side is not reached, we continue to explore in this direction,
+            // with a variable quality.
+            if(!reached1){
+                uv1 -= offset * QUALITY(i);
+            }
+            if(!reached2){
+                uv2 += offset * QUALITY(i);
+            }
+
+            // If both sides have been reached, stop the exploration.
+            if(reachedBoth){
+                break;
+            }
+        }
+    }
+
+    // Compute the distances to each side edge of the edge (!).
+    float distance1 = isHorizontal ? (coords.x - uv1.x) : (coords.y - uv1.y);
+    float distance2 = isHorizontal ? (uv2.x - coords.x) : (uv2.y - coords.y);
+
+    // In which direction is the side of the edge closer ?
+    bool isDirection1 = distance1 < distance2;
+    float distanceFinal = min(distance1, distance2);
+
+    // Thickness of the edge.
+    float edgeThickness = (distance1 + distance2);
+
+    // Is the luma at center smaller than the local average ?
+    bool isLumaCenterSmaller = lumaCenter < lumaLocalAverage;
+
+    // If the luma at center is smaller than at its neighbour,
+    // the delta luma at each end should be positive (same variation).
+    bool correctVariation1 = (lumaEnd1 < 0.0) != isLumaCenterSmaller;
+    bool correctVariation2 = (lumaEnd2 < 0.0) != isLumaCenterSmaller;
+
+    // Only keep the result in the direction of the closer side of the edge.
+    bool correctVariation = isDirection1 ? correctVariation1 : correctVariation2;
+
+    // UV offset: read in the direction of the closest side of the edge.
+    float pixelOffset = - distanceFinal / edgeThickness + 0.5;
+
+    // If the luma variation is incorrect, do not offset.
+    float finalOffset = correctVariation ? pixelOffset : 0.0;
+
+    // Sub-pixel shifting
+    // Full weighted average of the luma over the 3x3 neighborhood.
+    float lumaAverage = (1.0 / 12.0) * (2.0 * (lumaDownUp + lumaLeftRight) + lumaLeftCorners + lumaRightCorners);
+    // Ratio of the delta between the global average and the center luma,
+    // over the luma range in the 3x3 neighborhood.
+    float subPixelOffset1 = clamp(abs(lumaAverage - lumaCenter) / lumaRange, 0.0, 1.0);
+    float subPixelOffset2 = (-2.0 * subPixelOffset1 + 3.0) * subPixelOffset1 * subPixelOffset1;
+    // Compute a sub-pixel offset based on this delta.
+    float subPixelOffsetFinal = subPixelOffset2 * subPixelOffset2 * float(dSubpixelQuality);
+
+    // Pick the biggest of the two offsets.
+    finalOffset = max(finalOffset, subPixelOffsetFinal);
+
+    // Compute the final UV coordinates.
+    vec2 finalUv = coords;
+    if(isHorizontal){
+        finalUv.y += finalOffset * stepLength;
+    } else {
+        finalUv.x += finalOffset * stepLength;
+    }
+
+    // Read the color at the new UV coordinates, and use it.
+    gl_FragColor = texture2D(tColor, finalUv);
+}
+`;

src/mol-gl/webgl/resources.ts

@@ -22,7 +22,7 @@ import { VertexArray, createVertexArray } from './vertex-array';
 
 function defineValueHash(v: boolean | number | string): number {
     return typeof v === 'boolean' ? (v ? 1 : 0) :
-        typeof v === 'number' ? v : hashString(v);
+        typeof v === 'number' ? (v * 10000) : hashString(v);
 }
 
 function wrapCached<T extends Resource>(resourceItem: ReferenceItem<T>) {

src/mol-plugin-state/builder/structure/hierarchy-preset.ts

@@ -76,6 +76,7 @@ const defaultPreset = TrajectoryHierarchyPresetProvider({
 
 const AllModelsParams = (a: PluginStateObject.Molecule.Trajectory | undefined, plugin: PluginContext) => ({
     useDefaultIfSingleModel: PD.Optional(PD.Boolean(false)),
+    representationPresetParams: PD.Optional(PD.Group(StructureRepresentationPresetProvider.CommonParams)),
     ...CommonParams(a, plugin)
 });
 

src/mol-plugin/config.ts

@@ -24,7 +24,7 @@ export const PluginConfig = {
     General: {
         IsBusyTimeoutMs: item('plugin-config.is-busy-timeout', 750),
         DisableAntialiasing: item('plugin-config.disable-antialiasing', false),
-        PixelScale: item<number | undefined>('plugin-config.pixel-scale', void 0),
+        PixelScale: item('plugin-config.pixel-scale', 1),
         EnableWboit: item('plugin-config.enable-wboit', false)
     },
     State: {

src/mol-plugin/util/viewport-screenshot.ts

@@ -115,7 +115,10 @@ class ViewportScreenshotHelper extends PluginComponent {
                 mode: mutlisample ? 'on' : 'off',
                 sampleLevel: colorBufferFloat && textureFloat ? 4 : 2
             },
-            postprocessing: c.props.postprocessing
+            postprocessing: {
+                ...c.props.postprocessing,
+                antialiasing: { name: 'off', params: {} }
+            }
         });
     }
 
@@ -131,7 +134,10 @@ class ViewportScreenshotHelper extends PluginComponent {
                 cameraHelper: { axes: this.values.axes },
                 transparentBackground: this.values.transparent,
                 // TODO: optimize because this creates a copy of a large object!
-                postprocessing: this.plugin.canvas3d!.props.postprocessing
+                postprocessing: {
+                    ...this.plugin.canvas3d!.props.postprocessing,
+                    antialiasing: { name: 'off', params: {} }
+                }
             });
             return this._imagePass;
         }

src/mol-repr/structure/visual/util/link.ts

@@ -131,7 +131,7 @@ export function createLinkCylinderMesh(ctx: VisualContext, linkBuilder: LinkBuil
 
             addFixedCountDashedCylinder(builderState, va, vb, 0.5, 7, cylinderProps);
         } else if (linkStyle === LinkStyle.Double || linkStyle === LinkStyle.Triple) {
-            const order = LinkStyle.Double ? 2 : 3;
+            const order = linkStyle === LinkStyle.Double ? 2 : 3;
             const multiRadius = linkRadius * (linkScale / (0.5 * order));
             const absOffset = (linkRadius - multiRadius) * linkSpacing;