Skip to content
Snippets Groups Projects
Commit de8f2943 authored by David Sehnal's avatar David Sehnal
Browse files

alpha-orbitals extension and example

parent 013ddb72
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 60875 additions and 2 deletions
src/examples/alpha-orbitals/example-data.ts 0 → 100644
+ 60186
− 0
View file @ de8f2943
Source diff could not be displayed: it is too large. Options to address this: view the blob.
src/examples/alpha-orbitals/index.html 0 → 100644
+ 32
− 0
View file @ de8f2943
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, user-scalable=no, minimum-scale=1.0, maximum-scale=1.0">
<title>Mol* Alpha Orbitals Example</title>
<style>
* {
margin: 0;
padding: 0;
box-sizing: border-box;
}
#app {
position: absolute;
left: 160px;
top: 100px;
width: 600px;
height: 600px;
border: 1px solid #ccc;
}
</style>
<link rel="stylesheet" type="text/css" href="molstar.css" />
<script type="text/javascript" src="./index.js"></script>
</head>
<body>
<div id='controls'></div>
<div id="app"></div>
<script>
AlphaOrbitalsExample.init('app')
</script>
</body>
</html>
\ No newline at end of file
src/examples/alpha-orbitals/index.ts 0 → 100644
+ 99
− 0
View file @ de8f2943
/**
* Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
*
* @author Alexander Rose <alexander.rose@weirdbyte.de>
*/
import { Basis, computeIsocontourValues } from '../../extensions/alpha-orbitals/cubes';
import { SphericalBasisOrder } from '../../extensions/alpha-orbitals/orbitals';
import { createPluginAsync, DefaultPluginSpec } from '../../mol-plugin';
import { createVolumeRepresentationParams } from '../../mol-plugin-state/helpers/volume-representation-params';
import { StateTransforms } from '../../mol-plugin-state/transforms';
import { PluginContext } from '../../mol-plugin/context';
import { ColorNames } from '../../mol-util/color/names';
import { DemoMoleculeSDF, DemoOrbitals } from './example-data';
import './index.html';
import { CreateOrbitalVolume, StaticBasisAndOrbitals } from './transforms';
require('mol-plugin-ui/skin/light.scss');
interface DemoInput {
moleculeSdf: string,
basis: Basis,
order: SphericalBasisOrder,
orbitals: {
energy: number,
alpha: number[]
}[]
}
class AlphaOrbitalsExample {
plugin: PluginContext;
async init(target: string | HTMLElement) {
this.plugin = await createPluginAsync(typeof target === 'string' ? document.getElementById(target)! : target, {
...DefaultPluginSpec,
layout: {
initial: {
isExpanded: false,
showControls: false
},
controls: { left: 'none', right: 'none', top: 'none', bottom: 'none' }
}
});
this.plugin.managers.interactivity.setProps({ granularity: 'element' });
this.load({
moleculeSdf: DemoMoleculeSDF,
...DemoOrbitals
});
}
async load(input: DemoInput) {
await this.plugin.clear();
const data = await this.plugin.builders.data.rawData({ data: input.moleculeSdf }, { state: { isGhost: true } });
const trajectory = await this.plugin.builders.structure.parseTrajectory(data, 'mol');
const model = await this.plugin.builders.structure.createModel(trajectory);
const structure = await this.plugin.builders.structure.createStructure(model);
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 volumeRef = await this.plugin.build().toRoot()
.apply(StaticBasisAndOrbitals, { basis: input.basis, order: input.order, orbitals: input.orbitals })
.apply(CreateOrbitalVolume, { index: 44 })
.commit();
if (!volumeRef.isOk) return;
// TODO: the isovalues are being computed twice. Need to add more flexible support to Volume object
// for controlling them
const { negative, positive } = computeIsocontourValues(volumeRef.data!.grid.cells.data as any, 0.85);
const repr = this.plugin.build().to(volumeRef);
if (positive !== void 0) {
repr.apply(StateTransforms.Representation.VolumeRepresentation3D, createVolumeRepresentationParams(this.plugin, volumeRef.data!, {
type: 'isosurface',
typeParams: { isoValue: { kind: 'absolute', absoluteValue: positive } },
color: 'uniform',
colorParams: { value: ColorNames.blue }
}));
}
if (negative !== void 0) {
repr.apply(StateTransforms.Representation.VolumeRepresentation3D, createVolumeRepresentationParams(this.plugin, volumeRef.data!, {
type: 'isosurface',
typeParams: { isoValue: { kind: 'absolute', absoluteValue: negative } },
color: 'uniform',
colorParams: { value: ColorNames.red }
}));
}
await repr.commit();
}
}
(window as any).AlphaOrbitalsExample = new AlphaOrbitalsExample();
\ No newline at end of file
src/examples/alpha-orbitals/transforms.ts 0 → 100644
+ 74
− 0
View file @ de8f2943
/**
* Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
*
* @author David Sehnal <david.sehnal@gmail.com>
*/
import { PluginStateObject, PluginStateTransform } from '../../mol-plugin-state/objects';
import { Basis, createSphericalCollocationGrid } from '../../extensions/alpha-orbitals/cubes';
import { ParamDefinition as PD } from '../../mol-util/param-definition';
import { Task } from '../../mol-task';
import { CustomProperties } from '../../mol-model/custom-property';
import { SphericalBasisOrder } from '../../extensions/alpha-orbitals/orbitals';
import { Volume } from '../../mol-model/volume';
export class BasisAndOrbitals extends PluginStateObject.Create<{ basis: Basis, order: SphericalBasisOrder, orbitals: { energy: number, alpha: number[] }[] }>({ name: 'Basis', typeClass: 'Object' }) { }
export const StaticBasisAndOrbitals = PluginStateTransform.BuiltIn({
name: 'static-basis-and-orbitals',
display: 'Basis and Orbitals',
from: PluginStateObject.Root,
to: BasisAndOrbitals,
params: {
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 })
},
})({
apply({ params }) {
return new BasisAndOrbitals({ basis: params.basis, order: params.order, orbitals: params.orbitals }, { label: params.label });
}
});
export const CreateOrbitalVolume = PluginStateTransform.BuiltIn({
name: 'create-orbital-volume',
display: 'Orbital Volume',
from: BasisAndOrbitals,
to: PluginStateObject.Volume.Data,
params: (a) => {
if (!a) {
return { index: PD.Numeric(0) };
}
return {
index: PD.Select(0, a.data.orbitals.map((o, i) => [i, `[${i + 1}] ${o.energy.toFixed(4)}`]))
};
}
})({
apply({ a, params }) {
return Task.create('Orbital Volume', async ctx => {
const data = await createSphericalCollocationGrid({
basis: a.data.basis,
cutoffThreshold: 0.0075,
alphaOrbitals: a.data.orbitals[params.index].alpha,
sphericalOrder: a.data.order,
boxExpand: 4.5,
gridSpacing: [
[55, 0.6],
[40, 0.5],
[25, 0.4],
[0, 0.33],
],
}).runInContext(ctx);
const volume: Volume = {
grid: data.grid,
sourceData: { name: 'custom grid', kind: 'custom', data },
customProperties: new CustomProperties(),
_propertyData: Object.create(null),
};
return new PluginStateObject.Volume.Data(volume, { label: 'Orbital Volume' });
});
}
});
\ No newline at end of file
src/extensions/alpha-orbitals/collocation.ts 0 → 100644
+ 187
− 0
View file @ de8f2943
/**
* Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
*
* @author David Sehnal <david.sehnal@gmail.com>
*/
import { Vec3 } from '../../mol-math/linear-algebra';
import { RuntimeContext } from '../../mol-task';
import { arrayMin } from '../../mol-util/array';
import { Basis, CubeGridInfo, SphericalElectronShell } from './cubes';
import {
normalizeBasicOrder,
SphericalFunctions,
SphericalBasisOrder,
} from './orbitals';
export interface CollocationParams {
grid: CubeGridInfo;
basis: Basis;
sphericalOrder: SphericalBasisOrder;
cutoffThreshold: number;
alphaOrbitals: number[];
}
export async function sphericalCollocation(
{
grid,
basis,
sphericalOrder,
cutoffThreshold,
alphaOrbitals,
}: CollocationParams,
taskCtx: RuntimeContext
) {
let baseCount = 0;
for (const atom of basis.atoms) {
for (const shell of atom.shells) {
baseCount += 2 * shell.angularMomentum + 1;
}
}
const matrix = new Float32Array(grid.npoints);
let baseIndex = 0;
for (const atom of basis.atoms) {
for (const shell of atom.shells) {
const L = shell.angularMomentum;
if (L > 4) {
// TODO: will L > 4 be required? Would need to precompute more functions in that case.
throw new Error('Angular momentum L > 4 not supported.');
}
const alpha = normalizeBasicOrder(
L,
alphaOrbitals.slice(baseIndex, baseIndex + 2 * L + 1),
sphericalOrder
);
baseIndex += 2 * L + 1;
collocationBasis(
matrix,
grid,
atom.center,
cutoffThreshold,
alpha,
shell
);
if (taskCtx.shouldUpdate) {
await taskCtx.update({
message: 'Computing...',
current: baseIndex,
max: baseCount,
isIndeterminate: false,
});
}
}
}
return matrix;
}
function collocationBasis(
matrix: Float32Array,
grid: CubeGridInfo,
center: Vec3,
cutoffThreshold: number,
alpha: number[],
shell: SphericalElectronShell
) {
const L = shell.angularMomentum;
const exponents = shell.exponents;
const ncoeff = exponents.length;
const coefficients: number[] = sumCoefficients(shell.coefficients);
const sphericalFunc = SphericalFunctions[L];
const cx = center[0],
cy = center[1],
cz = center[2];
const ny = grid.dimensions[1],
nz = grid.dimensions[2];
const gdx = grid.delta[0],
gdy = grid.delta[1],
gdz = grid.delta[2];
const sx = grid.box.min[0],
sy = grid.box.min[1],
sz = grid.box.min[2];
const cutoffRadius =
cutoffThreshold > 0
? Math.sqrt(-Math.log(cutoffThreshold) / arrayMin(exponents))
: 10000;
const cutoffSquared = cutoffRadius * cutoffRadius;
const radiusBox = getRadiusBox(grid, center, cutoffRadius);
const iMin = radiusBox[0][0],
jMin = radiusBox[0][1],
kMin = radiusBox[0][2];
const iMax = radiusBox[1][0],
jMax = radiusBox[1][1],
kMax = radiusBox[1][2];
for (let i = iMin; i <= iMax; i++) {
const x = sx + gdx * i - cx;
const oX = i * ny * nz;
for (let j = jMin; j <= jMax; j++) {
const y = sy + gdy * j - cy;
const oY = oX + j * nz;
for (let k = kMin; k <= kMax; k++) {
const z = sz + gdz * k - cz;
const R2 = x * x + y * y + z * z;
if (R2 > cutoffSquared) {
continue;
}
let gaussianSum = 0;
for (let c = 0; c < ncoeff; c++) {
gaussianSum +=
coefficients[c] * Math.exp(-exponents[c] * R2);
}
const sphericalSum =
L === 0 ? alpha[0] : sphericalFunc(alpha, x, y, z);
matrix[k + oY] += gaussianSum * sphericalSum;
}
}
}
}
function sumCoefficients(coefficients: number[][]) {
if (coefficients.length === 1) return coefficients[0];
const ret: number[] = [];
const len = coefficients[0].length;
for (let i = 0; i < len; i++) ret[i] = 0;
for (let j = 0; j < coefficients.length; j++) {
const row = coefficients[j];
for (let i = 0; i < len; i++) ret[i] += row[i];
}
return ret;
}
function getRadiusBox(grid: CubeGridInfo, center: Vec3, radius: number) {
const r = Vec3.create(radius, radius, radius);
const min = Vec3.scaleAndAdd(Vec3(), center, r, -1);
const max = Vec3.add(Vec3(), center, r);
Vec3.sub(min, min, grid.box.min);
Vec3.sub(max, max, grid.box.min);
Vec3.div(min, min, grid.delta);
Vec3.floor(min, min);
Vec3.max(min, min, Vec3());
Vec3.div(max, max, grid.delta);
Vec3.ceil(max, max);
Vec3.min(max, max, Vec3.subScalar(Vec3(), grid.dimensions, 1));
return [min, max];
}
src/extensions/alpha-orbitals/cubes.ts 0 → 100644
+ 204
− 0
View file @ de8f2943
/**
* 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 { 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 { sphericalCollocation } from './collocation';
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[];
}[];
}
export interface SphericalElectronShell {
angularMomentum: number;
exponents: number[];
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[];
}
export function createSphericalCollocationGrid(
params: SphericalCollocationParams
): Task<CubeGrid> {
return Task.create('Spherical Collocation Grid', async (ctx) => {
const grid = initBox(
params.basis.atoms.map((a) => a.center),
params.gridSpacing,
params.boxExpand
);
const matrix = await sphericalCollocation(
{
grid,
basis: params.basis,
alphaOrbitals: params.alphaOrbitals,
cutoffThreshold: params.cutoffThreshold,
sphericalOrder: 'cca-reverse', // renamed entos ordering
},
ctx
);
return createCubeGrid(grid, matrix);
});
}
const BohrToAngstromFactor = 0.529177210859;
function createCubeGrid(gridInfo: CubeGridInfo, values: Float32Array) {
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, [0, 1, 2], Float32Array),
(values as any) as Tensor.Data
),
stats: {
min: arrayMin(values),
max: arrayMax(values),
mean: arrayMax(values),
sigma: arrayRms(values),
},
};
const 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(
values: Float32Array,
cumulativeThreshold: number
) {
const size = values.length;
const weights = new Float32Array(size);
const indices = new Int32Array(size);
let weightSum = 0;
for (let i = 0; i < size; i++) {
const v = values[i];
const w = v * v;
weights[i] = w;
indices[i] = i;
weightSum += w;
}
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/orbitals.ts 0 → 100644
+ 91
− 0
View file @ de8f2943
/**
* Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
*
* @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)
);
}
webpack.config.js
+ 1
− 1
View file @ de8f2943
const { createApp, createExample, createBrowserTest } = require('./webpack.config.common.js'); const { createApp, createExample, createBrowserTest } = require('./webpack.config.common.js');
const examples = ['proteopedia-wrapper', 'basic-wrapper', 'lighting']; const examples = ['proteopedia-wrapper', 'basic-wrapper', 'lighting', 'alpha-orbitals'];
const tests = [ const tests = [
'font-atlas', 'font-atlas',
'marching-cubes', 'marching-cubes',
......
webpack.config.production.js
+ 1
− 1
View file @ de8f2943
const { createApp, createExample } = require('./webpack.config.common.js'); const { createApp, createExample } = require('./webpack.config.common.js');
const examples = ['proteopedia-wrapper', 'basic-wrapper', 'lighting']; const examples = ['proteopedia-wrapper', 'basic-wrapper', 'lighting', 'alpha-orbitals'];
module.exports = [ module.exports = [
createApp('viewer', 'molstar'), createApp('viewer', 'molstar'),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment