Изгиб плоской молекулы

[B@R5uk]

Тогда проблема.

С другой стороны, поскольку смещение по Z мало, то поправки на координаты X и Y будут вообще второго порядка малости. Ими точно нельзя пренебречь?

-- 06.08.2019, 23:37 --

madschumacher, у вас есть возможность с нуля рассчитать координаты атомов в плоскости? Или они берутся "из справочника" в каждом случае?

[madschumacher]

B@R5uk в сообщении #1409012 писал(а):

Тогда проблема.

Я заметил. Поэтому и разместил вопрос в ПРР(М).

B@R5uk в сообщении #1409012 писал(а):

Ими точно нельзя пренебречь?

Ну наверное, раз я пишу, что нельзя, наверное я уже пробовал, и знаю, что нельзя?

B@R5uk в сообщении #1409012 писал(а):

у вас есть возможность с нуля рассчитать координаты атомов в плоскости?

Да, могу и рассчитать, могу и из справочника (базы данных) стырить.

[B@R5uk]

madschumacher в сообщении #1409013 писал(а):

Да, могу и рассчитать...

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

[madschumacher]

B@R5uk в сообщении #1409016 писал(а):

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

Ага, спасибо. Если бы это было так просто, я б не спрашивал об извращенных аналитичесих методах отображения плоскости на более сложную поверхность. Тем более, про простую версию этого, предложенную DeBill (модель пружинок) я уже комментировал, говоря, что это существенно усложняет задачу.

B@R5uk в сообщении #1409016 писал(а):

Длины связей должны остаться теми же, а вот углы между ними поплывут.

А это уже другое решение. У Вас есть идеи как сделать это конкретно, алгоритмически?

[B@R5uk]

madschumacher в сообщении #1409006 писал(а):

Об этом я тоже писал, и сказал, что некие деформации все ж допустимы...

Я так понимаю, вы имеете в виду случаи цилиндрической деформации. Этому есть простое объяснение: у плоскости и цилиндра одна и та же кривизна, поэтому плоскость можно навернуть на цилиндр без растяжений. В общем случае этого не получится, поэтому надо честно считать.

-- 06.08.2019, 23:56 --

madschumacher в сообщении #1409017 писал(а):

У Вас есть идеи как сделать это конкретно, алгоритмически?

А вы можете поделиться тем, что у вас уже считает в плоском случае? Или секрет фирмы? Просто не хочется строить алгоритм на пустом месте.

[madschumacher]

B@R5uk в сообщении #1409018 писал(а):

А вы можете поделиться тем, что у вас уже считает в плоском случае? Или секрет фирмы? Просто не хочется строить алгоритм на пустом месте.

import numpy as np
from math import *
from copy import deepcopy

inpf=open("molecule.xyz", "r")

atn=[]
xyz=[]

nline=0
nat=0
for line in inpf:
    nline+=1
    words=line.split()
    if nline==1:
        nat=int(words[0])
    elif nline>2 and len(words)>3:
        atn.append(words[0])
        xyz.append(np.array(words[1:4], dtype=float))

if nat!=len(xyz):
    print("Number of atoms specified is inconsistent with xyz definition")
    exit(1)

xix=np.array([0.0,  0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2, 0.0, -0.2, -0.4, -0.6, -0.8, -1.0])*0.05
xiy=np.array([0.0 for i in xix])
def new_xyz(r, txix, txiy):
    tr=deepcopy(r)
    tr[-1] = txix*tr[0]**2 + txiy*tr[1]**2

    nold=np.array([0.0, 0.0, -1.])
    nold/=sqrt(np.dot(nold,nold))
    nnew=np.array([2.0*txix*tr[0], 2.0*txiy*tr[1], -1.0])
    nnew/=sqrt(np.dot(nnew,nnew))

    rpr=np.dot(nold,r)
    
    return (tr+nnew*rpr)

for i in range(0,len(xix)):
    print("%i\n %f  %f" % (nat, xix[i], xiy[i]))
    for n in range(0,len(atn)):
        tr=new_xyz(xyz[n], xix[i], xiy[i])
        print(" %3s  %15.4f %15.4f %15.4f "  % (atn[n],  tr[0], tr[1], tr[2]))

[ИСН]

Не надо извращённых методов. На поверхность, которую можно свернуть из плоскости, плоскость отображается легко. А на такую, которую нельзя - не отображается ни легко, ни сложно. Иначе зачем бы, думаете, люди морочились с проекцией Меркатора и другими? Взяли бы да развернули глобус на плоскость. Но нет, от искажений так не уйти.

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

[пианист]

madschumacher в сообщении #1408993 писал(а):

Да, очень похоже на то, что хочется. Ну или с изгибами листа бумаги (с изломами, которые будут убираться, как на мыльной плёнке).
Это можно как-то описать в человеческих терминах?

Как свойства решений задачи Плато.
Кстати, описание асимптотик по малым отклонениям контура может оказаться довольно забавным, независимо от.

[madschumacher]

ИСН в сообщении #1409023 писал(а):

Иначе зачем бы, думаете, люди морочились с проекцией Меркатора и другими?

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

ИСН в сообщении #1409023 писал(а):

Но нет, от искажений так не уйти.

Ну я об этом и говорю: искажение - - это не страшно, лишь бы они были меньше, чем в простом отображении на поверхность 2го порядка.

ИСН в сообщении #1409023 писал(а):

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

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

пианист в сообщении #1409031 писал(а):

Как свойства решений задачи Плато.

Это та, о которой Фоменко как-то в Кванте писал? А где можно почитать про это, чтобы язык был понятен не-математикам?

[B@R5uk]

madschumacher в сообщении #1409021 писал(а):

Код:

Я имел в виду не текст кода, а именно алгоритм словами.

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

ИСН в сообщении #1409023 писал(а):

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

По хорошему то, да... вот только это сложная задача. И подход топикстартера может как раз к этому и ведёт: если в базисе малых искажений рассчитать коэффициенты жёсткости (в том числе перекрёстные), то можно брать различные начальные условия (деформации) и рассчитывать динамику того, как они друг в друга переходят. Потом посмотреть на спектры колебаний каждого коэффициента. Всё, что соответствует одной частоте — будет нормальное колебание. Таким образом нормальное колебание будет разложено по базису простых малых деформаций. Но тут самое главное, чтобы для частот интереса базис оказался полным.

[madschumacher]

B@R5uk в сообщении #1409039 писал(а):

Я имел в виду не текст кода, а именно алгоритм словами.

madschumacher в сообщении #1408958 писал(а):

Банальное отображение на плоскость второго порядка, так что координаты $(x,y,0) \rightarrow (x,y,z)$, где $z=k_x x^2 + k_y y^2 + k_{xy} xy$ не канают

[пианист]

В Кванте не читал, но у Фоменко, да, есть результаты по этой теме, так что вполне может быть.
Я лично читал про Плато по несколько странному источнику ;), есть такая книжечка (в сети имеется), Полищук "Софус Ли", вот там.
Гляньте, вроде все понятно изложено.

[B@R5uk]

Ладно. Мне видится только такое решение задачи. На координаты всех атомов накладывается ограничение, чтобы они лежали на поверхности 2-го порядка. На длины всех связей накладывается ограничение, чтобы они оставались неизменными. Из всех углов между связями вычитаются их значения в плоской молекуле, разности возводятся в квадрат и суммируются с некими весами, характеризующими жёсткость соответствующих углов между связями. В результате получается некая целевая функция и набор ограничений. Численно рассчитываем минимум этой целевой функции при заданных ограничениях. Получаем искомые координаты.

Это общий план, в него можно внести множество оптимизаций. Например, координаты $\{x,y,z\}$ будут весьма неудачными для сложности расчёта, их лучше заменить на углы, поскольку длины связей фиксированы. Ограничения на длину связи лучше не использовать в расчётах в лоб, а задаться небольшим числом новых угловых координат и по топологии молекулы рассчитать новые значения координат атомов. Это сразу удовлетворит требуемым ограничениям, и о них в последствии не нужно будет беспокоиться. Затем это всё подставляется в целевую функцию и оптимизируется по новым координатам.

[upgrade]

Эту задачу можно так сформулировать?:
$s=\int\limits_{a}^{b} \sqrt{x'^2(t)+y'^2(t)+z'^2(t)}dt$
Известно $s$ (длина расстояния между точками) и вид функции, заданной параметрически (форма изгиба связей между точками).
Необходимо найти $x,y,z$

[madschumacher]

upgrade в сообщении #1409049 писал(а):

Эту задачу можно так сформулировать?:

Да, мне как раз тоже пришла такая идея.

В итоге реализовал что-то, что вроде меня сейчас устраивает. Привожу то, что сделал для проверки.
Допустим, у нас есть изначальный вектор $\mathbf{r}_0 = (x_0, y_0, 0)^\dagger$, мы его хотим перевести в вектор, который параметризован как $\mathbf{r} = (x_0 t, y_0 t , z(t) )^\dagger$, где $z = \underbrace{(ax_0^2 + bx_0 y_0 + cy_0^2)}_{q} t^2$ (т.е. мы уменьшаем проекцию $\mathbf{r}$ на плоскость $xy$, но эта проекция остаётся в том же направлении).
Дополнительным условием сделаем, что расстояние от начала координат $(0,0,0)^\dagger$ до новой точки по поверхности должно остаться прежним, т.е.
$\int \limits_0^{t_\text{fin}} \sqrt{\underbrace{x_0^2 + y_0^2}_{A} + \underbrace{(2q)^2}_{B}t^2} = \frac{A \arcsinh\left(\sqrt{\frac{B}{A}} t_\text{fin}\right)}{2 \sqrt{B}} + \frac{ t_\text{fin} \sqrt{B t_\text{fin}^2 + A}}{2} = |\mathbf{r}_0| = \sqrt{x_0^2 + y_0^2}$.
Решая это уравнение, находим $t_\text{fin}$, которое и даёт нам искомые координаты.
На примере коронена получилось вот что:

(код этого алгоритма)

import numpy as np
from math import *
from copy import deepcopy
import scipy.optimize as spo

inpf=open("coronene.xyz", "r")

atn=[]
xyz=[]

nline=0
nat=0
for line in inpf:
    nline+=1
    words=line.split()
    if nline==1:
        nat=int(words[0])
    elif nline>2 and len(words)>3:
        atn.append(words[0])
        xyz.append(np.array(words[1:4], dtype=float))

if nat!=len(xyz):
    print("Number of atoms specified is inconsistent with xyz definition")
    exit(1)


def curv_length(t, A, B):
    return (A*np.arcsinh(sqrt(B/A)*t)/(2.0*sqrt(B)) + 0.5*t*sqrt(B*t**2 + A) )

def lsqf(t, r0, A, B):
    return (curv_length(t, A, B) - r0)**2

def get_new_coords(xyz, a, b, c):
    lvec=sqrt(np.dot(xyz,xyz))
    newxyz = deepcopy(xyz)
    Q=2.0*(a*xyz[0]**2 + b*xyz[0]*xyz[1] + c*xyz[1]**2)
    if lvec>1.0e-2 and abs(Q)>1.0e-2:
        A=xyz[0]**2 + xyz[1]**2
        B=Q**2
        res=spo.minimize_scalar(lsqf, bounds=(0., 1.), args=(lvec, A,B), method='Golden')
        if res.success:
            newxyz[0]=xyz[0]*res.x
            newxyz[1]=xyz[1]*res.x
            newxyz[2]=0.5*Q*res.x**2
        else:
            print("Failed convergence in coordinates!")
            exit(1)
    return newxyz

for i in range(0,100):
    abc=0.05*2.0*(np.random.rand(3)-np.array([0.5,0.5,0.5]))
    print("%i\n %f  %f  %f" % (nat, abc[0], abc[1], abc[2]))
    for n in range(0,len(atn)):
        tr=get_new_coords(xyz[n], abc[0], abc[1], abc[2])
        print(" %3s  %15.4f %15.4f %15.4f "  % (atn[n],  tr[0], tr[1], tr[2]))
    

Спасибо большое всем за помощь и участие!