Как мы SciArt-ом логотип конференции прокачали

24.03.2021

В прошлом (2020) году в связи с пандемией мы проводили научную онлайн конференцию по вычислительной химии, и для неё сделали логотип, который был, мягко говоря, так себе. Под катом рассказ о том, как мы его прокачали для конференции этого (2021) года при помощи небольшого количества квантовой механики, метода Монте-Карло, Python и Gnuplot.

Немного предыстории

2020 год, коронавирус набирает силу, вводятся ограничения, и для поддержания уровня научной дискуссии, улучшения научного сотрудничества среди русскоязычных учёных, разбросанных по всему миру, и для всего остального хорошего, мы организовывали в онлайне мини-симпозиум по вычислительной химии.

Но что это за конференция без логотипа? И собрав в кулак все художественные (от слова «худо») возможности оргкомитета конференции, в Inkscape был выстрадан логотип:

Как мы SciArt-ом логотип конференции прокачали — Изначальный логотип, который мы сделали при помощи Inkscape и Jmol+Gimp (для рисования фуллеренов).

Содержание картинки простое. В верхней левой части у нас есть светлячок, который по совместительству d-образная орбиталь. Это референс к квантово-химическому пакету Firefly, его создатель, Александр Александрович Грановский, безвременно ушёл из жизни за несколько месяцев до этого, и конференция была посвящена ему. Справа летают бакминстерфуллерены, слева внизу присутствует стилизованный гауссовообразный лазерный импульс. Ну и в самом низу олицетворение всея квантовой механики: буква ψ, привычно обозначающая волновую функцию, эта буква стилизована под факел. Несмотря на присутствие огня внутри картинки, сам логотип получился далеко не огонь. Поэтому в этом году мы предприняли попытку to Pimp My ~~Ride~~ Logo.

Генерируем светлячка из настоящих орбиталей

Основным элемента логотипа, который был в 2020 году, был светлячок-орбиталь. Его мы и захотели сохранить и улучшить. И что может быть более крутым, чем настоящая квантово-механическая орбиталь?

Поэтому мы решили собрать нового светлячка из настоящих орбиталей. Ограничиться решили несколькими составными частями:

крылья,
туловище,
голова,
и последним элементом было свечение туловища, огонь.

Все эти штуки мы решили представить в виде каких-то орбиталей ψ (волновых функций), зависящих от координат (x,y,z)=r. Красивой и простейшей идеей показалось отобразить квантовость при помощи семплирования по методу Метрополиса. Волновая функция в данном случае может быть представлена как траектория точек r₀, r₁, r₂, …, r_N, N — это длина траектории. В каждой из этих точек нам ещё известно значение волновой функции ψ_n=ψ(r_n). Сам алгоритм семплирования простейший, и пишется на Python за три минуты:

находясь в точке r_n мы генерируем новую пробную точку r_trial, добавляя к каждой координате (x_n,y_n,z_n) точки r_n случайную добавку в диапазоне от -δ до +δ, где δ — некое число, которое мы подбирали для каждого из случаев отдельно,
вероятность принятия новой точки равна p_acc = |ψ_trial|²/|ψ_n|², генерируя значение q из равномерного распределения от 0 до 1, мы сравниваем p_acc и q: если p_acc<q, то мы остаёмся в той же точке (r_n₊₁=r_n), а если нет, то переходим в новую точку (r_n₊₁=r_trial),
в случае если мы остались в той же точке, чтобы не было скучно и было больше точек, мы добавляем к сохраняемой в траектории точке случайное смещение по каждой координате в диапазаоне от -0.1⋅δ до +0.1⋅δ.

Стартовать траекторию, естественно, надо с точки с ненулевой волновой функцией. Первые же 10% траектории мы выкидывали, но, поскольку симуляции быстрые, то это не особо вычислительно напряжно. Количество точек для каждого куска светлячка подбирали отдельно.

Крылья

Начать решили с крыльев. В качестве прообраза взяли d-орбиталь вида:

r — это модуль вектора r, нормировка для алгоритма Метрополиса не важна, поэтому здесь и далее мы её опускаем. Это выражение даёт четырёхлистник с одинаковыми лепестками, что не очень похоже на бабочку. Чтобы исправить это недоразумение, мы немного модифицировали эту функцию, добавив немного асимметрии:

Из необычного, мы сохраняли положительные и отрицательные значения волновой функции в разные файлы, чтобы строить их отдельно.

Код на Python для генерации крыльев

import numpy as np
import random as rnd

DXYZ = 4.0
dxyz = 0.1*DXYZ

def getWFN(R, a=0.05, b=0.9, R0=1.0):
    return (a*(R[1])**3 + b*R[0]*R[1])*np.exp(-np.sqrt(np.dot(R,R))/R0)

Npts = 40000
N2ignore = Npts/10

outf = open("wfn.dat", "w")
outfp = open("wfn_pos.dat", "w")
outfn = open("wfn_neg.dat", "w")

genNewXYZ = lambda xyz, shift: xyz + np.array([rnd.uniform(-shift,shift) for q in range(0,3)])

XYZ = np.array([1.0,1.0,0.0])
Psi = getWFN(XYZ, DXYZ)

for i in range(0,Npts):
    trialXYZ = genNewXYZ(XYZ, DXYZ)
    trialPsi = getWFN(trialXYZ)
    if i % N2ignore == 0:
        print("step #%i" % (i))

    Ptrial = rnd.random()
    Ptest = trialPsi**2/Psi**2

    Accepted = False
    if Ptrial < Ptest:
        Psi = trialPsi
        XYZ = trialXYZ
        Accepted = True

    if i>N2ignore:
        r2save = XYZ
        if not Accepted:
            r2save = genNewXYZ(XYZ, dxyz)
        outf.write(" %15.10f %15.10f %15.10f %15.5en" % tuple(list(r2save)+[Psi]))
        if Psi > 0.0:
            outfp.write(" %15.10f %15.10f %15.10f %15.5en" % tuple(list(r2save)+[Psi]))
        else:
            outfn.write(" %15.10f %15.10f %15.10f %15.5en" % tuple(list(r2save)+[Psi]))

Туловище

Для построения туловища была выбрана более сложная конструкция:

$psi(x,y,z) =begin{cases} (y^3 + 0.5cdot y)cdot expleft( -frac{|x|}{0.2}-frac{|y|}{0.5} -frac{|z|}{0.2}right) , y<0 \ 0 , ygeq 0 end{cases}$

Код на Python для генерации туловища

import numpy as np
import random as rnd

DXYZ = 0.5
dxyz = 0.1*DXYZ

def getWFN(R, a=1.0, b=0.5, X0=0.2, Y0=0.5, Z0=0.2):
    if R[1]<0.:
        return (a*(R[1])**3 + b*R[1])*np.exp(-np.abs(R[0])/X0 - np.abs(R[1])/Y0 - np.abs(R[2])/Z0)
    else: 
        return 0.0

Npts = 20000
N2ignore = Npts/10

outf = open("body.dat", "w")

genNewXYZ = lambda xyz, shift: xyz + np.array([rnd.uniform(-shift,shift) for q in range(0,3)])

XYZ = np.array([0.0,-0.2,0.0])
Psi = getWFN(XYZ, DXYZ)

for i in range(0,Npts):
    trialXYZ = genNewXYZ(XYZ, DXYZ)
    trialPsi = getWFN(trialXYZ)
    if i % N2ignore == 0:
        print("step #%i" % (i))

    Ptrial = rnd.random()
    Ptest = trialPsi**2/Psi**2

    Accepted = False
    if Ptrial < Ptest:
        Psi = trialPsi
        XYZ = trialXYZ
        Accepted = True

    if i>N2ignore:
        r2save = XYZ
        if not Accepted:
            r2save = genNewXYZ(XYZ, dxyz)

        outf.write(" %15.10f %15.10f %15.10f %15.5en" % tuple(list(r2save)+[Psi]))

Голова

Голова была представлена функцией, напоминающей s-орбиталь:

$psi(x,y,z) = expleft(-frac{|x|}{0.3}-frac{|y|}{0.5}-frac{|z|}{0.3} right)$

Код на Python для генерации головы

import numpy as np
import random as rnd

DXYZ = 0.9
dxyz = 0.1*DXYZ

def getWFN(R, a=1.0, b=0.5, X0=0.3, Y0=0.5, Z0=0.3):
    return np.exp(-np.abs(R[0])/X0 - np.abs(R[1])/Y0 - np.abs(R[2])/Z0)

Npts = 10000
N2ignore = Npts/10

outf = open("head.dat", "w")

genNewXYZ = lambda xyz, shift: xyz + np.array([rnd.uniform(-shift,shift) for q in range(0,3)])

XYZ = np.array([0.0,-0.2,0.0])
Psi = getWFN(XYZ, DXYZ)

for i in range(0,Npts):
    trialXYZ = genNewXYZ(XYZ, DXYZ)
    trialPsi = getWFN(trialXYZ)
    if i % N2ignore == 0:
        print("step #%i" % (i))

    Ptrial = rnd.random()
    Ptest = trialPsi**2/Psi**2

    Accepted = False
    if Ptrial < Ptest:
        Psi = trialPsi
        XYZ = trialXYZ
        Accepted = True

    if i>N2ignore:
        r2save = XYZ
        if not Accepted:
            r2save = genNewXYZ(XYZ, dxyz)

        outf.write(" %15.10f %15.10f %15.10f %15.5en" % tuple(list(r2save)+[Psi]))

Огонь вокруг туловища

Огонь был сгенерирован функцией того же вида, что и туловище, но со слегка подкрученными параметрами:

$psi(x,y,z) =begin{cases} (y^3 + 0.5cdot y)cdot expleft( -frac{|x|}{0.4}-frac{|y|}{0.68} -frac{|z|}{0.3}right) , y<0 \ 0 , ygeq 0 end{cases}$

Код на Python для генерации огня вокруг туловища

import numpy as np
import random as rnd

DXYZ = 0.9
dxyz = 0.1*DXYZ

def getWFN(R, a=1.0, b=0.5, X0=0.4, Y0=0.68, Z0=0.3):
    if R[1]<0.:
        return (a*(R[1])**3 + b*R[1])*np.exp(-np.abs(R[0])/X0 - np.abs(R[1])/Y0 - np.abs(R[2])/Z0)
    else: 
        return 0.0

Npts = 1000
N2ignore = Npts/10

outf = open("body_fire.dat", "w")

genNewXYZ = lambda xyz, shift: xyz + np.array([rnd.uniform(-shift,shift) for q in range(0,3)])

XYZ = np.array([0.0,-0.2,0.0])
Psi = getWFN(XYZ, DXYZ)

for i in range(0,Npts):
    trialXYZ = genNewXYZ(XYZ, DXYZ)
    trialPsi = getWFN(trialXYZ)
    if i % N2ignore == 0:
        print("step #%i" % (i))

    Ptrial = rnd.random()
    Ptest = trialPsi**2/Psi**2

    Accepted = False
    if Ptrial < Ptest:
        Psi = trialPsi
        XYZ = trialXYZ
        Accepted = True

    if i>N2ignore:
        r2save = XYZ
        if not Accepted:
            r2save = genNewXYZ(XYZ, dxyz)

        outf.write(" %15.10f %15.10f %15.10f %15.5en" % tuple(list(r2save)+[Psi]))

Рисуем светлячка

Питоновские скрипты сгенерировали нам метрополисовские траектории вида (x_n, y_n, z_n, ψ_n) для n=0,1,2,…, и их надо было превратить в картинку. Для этого мы использовали Gnuplot в режиме multiplot. В качестве вида мы взяли проекцию карты (set view map). Основной задачей оказался подобор палитр для раскрашивания бабочки, их мы выбирали из представленных в этом репозитории на Гитхабе: https://github.com/Gnuplotting/gnuplot-palettes.

Для тела очень быстро подобралась палитра inferno, для огонька — ylrd, а вот с крылями пришлось возиться больше. Из более-менее применимых было три варианта: plasma, parula, и итоговый вариант — prgn.

Plasma — классная тема, но фиолетовый на краях очень сильно сливался с тёмно-синим фоном, а огонёк туловища был слабо заметен в паре с жёлтым крылом. Parula была хороша всем, кроме того, что вызывала у некоторых членов оргкомитета ассоциации с цветами …шведского флага, поэтому, чтобы лишний раз не провоцировать политических флеймов вокруг научного мероприятия, от этой схемы мы отказались. В итоге удалось подобрать схему prgn, которая хороша всем: не сливается с фоном, огонёк хорошо заметен, да и сама схема приятна глазу. Ну а дальше манипуляциями в Inkscape получился окончательный вариант логотипа.

Скрипт для построения светлячка в Gnuplot

#set terminal pngcairo size 1500,1500   transparent 
set terminal pngcairo size 1500,1500  background rgb '#080045' 
set output "ff_v2.png" 

set multiplot

set view map 
set size ratio 1

unset colorbox
unset border
unset key
unset xtics
unset ytics

XMax = 7.

set xrange [-XMax:XMax]
set yrange [-XMax:XMax]

set pm3d depthorder hidden3d 1
set hidden3d

set style fill  transparent solid 0.35 noborder
set style circle radius 0.02

load 'prgn.pal'

#load 'plasma.pal'
#load 'parula.pal'
splot 'wfn_pos.dat' u 2:1:3:4 w p pt 7 ps 2 lc palette


splot 'wfn_neg.dat' u 2:1:3:4 w p pt 7 ps 2 lc palette
      
load 'inferno.pal'
splot 'body.dat' u 1:2:3:4 w l lc palette
      
splot 'head.dat' u 1:2:3:4 w l lc palette

load 'ylrd.pal'
splot 'body_fire.dat' u 1:2:3:(-$4)  w l lw 3 lc palette

Заключение

Рисование картинок из квантовой механики — забавная штука (см. например ещё этот пост): не очень простая, гемморойная, но результат обычно стоит затраченных усилий. 🙂

P.S.

Сомневаюсь, что кому-то будет особо интересно смотреть материалы конференции, но вдруг. Поэтому оставлю ссылку на её прошлогодний сайт, на сайт этого года, а также на YouTube канал с записями докладов.

Источник

Gnuplot, python, sciart, квантовая физика, монте-карло