В прошлом (2020) году в связи с пандемией мы проводили научную онлайн конференцию по вычислительной химии, и для неё сделали логотип, который был, мягко говоря, так себе. Под катом рассказ о том, как мы его прокачали для конференции этого (2021) года при помощи небольшого количества квантовой механики, метода Монте-Карло, Python и Gnuplot.
Немного предыстории
2020 год, коронавирус набирает силу, вводятся ограничения, и для поддержания уровня научной дискуссии, улучшения научного сотрудничества среди русскоязычных учёных, разбросанных по всему миру, и для всего остального хорошего, мы организовывали в онлайне мини-симпозиум по вычислительной химии.
Но что это за конференция без логотипа? И собрав в кулак все художественные (от слова "худо") возможности оргкомитета конференции, в Inkscape был выстрадан логотип:
Содержание картинки простое. В верхней левой части у нас есть светлячок, который по совместительству d-образная орбиталь. Это референс к квантово-химическому пакету Firefly, его создатель, Александр Александрович Грановский, безвременно ушёл из жизни за несколько месяцев до этого, и конференция была посвящена ему. Справа летают бакминстерфуллерены, слева внизу присутствует стилизованный гауссовообразный лазерный импульс. Ну и в самом низу олицетворение всея квантовой механики: буква ψ, привычно обозначающая волновую функцию, эта буква стилизована под факел. Несмотря на присутствие огня внутри картинки, сам логотип получился далеко не огонь. Поэтому в этом году мы предприняли попытку to Pimp My Ride Logo.
Генерируем светлячка из настоящих орбиталей
Основным элемента логотипа, который был в 2020 году, был светлячок-орбиталь. Его мы и захотели сохранить и улучшить. И что может быть более крутым, чем настоящая квантово-механическая орбиталь?
Поэтому мы решили собрать нового светлячка из настоящих орбиталей. Ограничиться решили несколькими составными частями:
крылья,
туловище,
голова,
и последним элементом было свечение туловища, огонь.
Все эти штуки мы решили представить в виде каких-то орбиталей ψ (волновых функций), зависящих от координат (x,y,z)=r. Красивой и простейшей идеей показалось отобразить квантовость при помощи семплирования по методу Метрополиса. Волновая функция в данном случае может быть представлена как траектория точек r0, r1, r2, ..., rN, N - это длина траектории. В каждой из этих точек нам ещё известно значение волновой функции ψn=ψ(rn). Сам алгоритм семплирования простейший, и пишется на Python за три минуты:
находясь в точке rn мы генерируем новую пробную точку rtrial, добавляя к каждой координате (xn,yn,zn) точки rn случайную добавку в диапазоне от -δ до +δ, где δ — некое число, которое мы подбирали для каждого из случаев отдельно,
вероятность принятия новой точки равна pacc = |ψtrial|2/|ψn|2, генерируя значение q из равномерного распределения от 0 до 1, мы сравниваем pacc и q: если pacc<q, то мы остаёмся в той же точке (rn+1=rn), а если нет, то переходим в новую точку (rn+1=rtrial),
в случае если мы остались в той же точке, чтобы не было скучно и было больше точек, мы добавляем к сохраняемой в траектории точке случайное смещение по каждой координате в диапазаоне от -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.5e\n" % tuple(list(r2save)+[Psi]))
if Psi > 0.0:
outfp.write(" %15.10f %15.10f %15.10f %15.5e\n" % tuple(list(r2save)+[Psi]))
else:
outfn.write(" %15.10f %15.10f %15.10f %15.5e\n" % tuple(list(r2save)+[Psi]))
Туловище
Для построения туловища была выбрана более сложная конструкция:
Код на 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.5e\n" % tuple(list(r2save)+[Psi]))
Голова
Голова была представлена функцией, напоминающей s-орбиталь:
Код на 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.5e\n" % tuple(list(r2save)+[Psi]))
Огонь вокруг туловища
Огонь был сгенерирован функцией того же вида, что и туловище, но со слегка подкрученными параметрами:
Код на 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.5e\n" % tuple(list(r2save)+[Psi]))
Рисуем светлячка
Питоновские скрипты сгенерировали нам метрополисовские траектории вида (xn, yn, zn, ψ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 канал с записями докладов.