расчет движения электронов в атоме лития методом МД(2мерная)

[MathKvant]

Привет всем, возник вопрос по МД. Задача состоит в нахождении траектории движения для каждого электрона в атоме лития. Для нахождения скорости и положения в пространстве использую центрально-разностную апроксимацию.
В итоге имею вот такие уравнения:
$v_x(t+2\tau)=v_x(t)+\frac{2 \tau F}{m} \\ X(t+2\tau)=X(t)+2\tau v_x(t+\tau) $
Аналогичные для y. $\tau$-шаг по времени, равный 0,01 от периода вращения вокруг ядра.
Как найти $v_x(t+\tau)$? Я делал типо такого $v_x(t+\tau)=(v_x(t+2\tau)+v_x(t))/2$, но не знаю верно ли? Вместо F подставляю сумму всех кулоновских сил спроецированные на оси x и y соответственно. Также не знаю какие x и y нужно подставлять для начала подсчета.

[Pphantom]

MathKvant в сообщении #972193 писал(а):

Привет всем, возник вопрос по МД. Задача состоит в нахождении траектории движения для каждого электрона в атоме лития.

Во-первых, что такое "МД"? Молекулярная динамика?

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

MathKvant в сообщении #972193 писал(а):

Для нахождения скорости и положения в пространстве использую центрально-разностную апроксимацию.

И как, получилось? Собственно, первый вопрос, который тут появляется - зачем Вы ее используете?

[MathKvant]

Pphantom в сообщении #972201 писал(а):

Во-первых, что такое "МД"? Молекулярная динамика?

Да, оно самое

Pphantom в сообщении #972201 писал(а):

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

Не я придумал задачу, может и глупо рассматривать траектория движения электрона, я это представляю как очень маленькие планеты.

Pphantom в сообщении #972201 писал(а):

И как, получилось? Собственно, первый вопрос, который тут появляется - зачем Вы ее используете?

$\\ \frac{dv}{dt}=\frac{F}{m}\\ \frac{dx}{dt}=v\\ f'(t)=\frac{f(t+\tau)-f(t-\tau)}{2\tau} $
Такое поставили задание, вот и использую

[Pphantom]

MathKvant в сообщении #972209 писал(а):

Не я придумал задачу, может и глупо рассматривать траектория движения электрона, я это представляю как очень маленькие планеты.

Ну хорошо, будем считать, что это данность.

MathKvant в сообщении #972209 писал(а):

\\$\frac{dv}{dt}=\frac{F}{m}\\ \frac{dx}{dt}=v\\ f'(t)=\frac{f(t+\tau)-f(t-\tau)}{2\tau} $

Вопрос не в том, что такое центрально-разностная аппроксимация (я это знаю), а в том, зачем Вы ее используете. Просто если Вы пытаетесь таким образом осложнить себе жизнь во имя какой-то великой цели , то это одно, а если это было сделано "просто так", то проще просто так не делать.

[MathKvant]

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

[Pphantom]

MathKvant в сообщении #972215 писал(а):

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

На форуме принято обращение на "Вы". И чем Вам односторонняя аппроксимация не нравится?

[мат-ламер]

MathKvant
Учебники по квантовой химии пробовали открывать? В лоб решать задачу будет затруднительно, поскольку волновая функция многомерная. Лучше всё-таки что сделано освоить. Допустим, симметрия сокращает вычисления. Но я в этом не спец.

[DimaM]

мат-ламер в сообщении #972255 писал(а):

Учебники по квантовой химии пробовали открывать? В лоб решать задачу будет затруднительно, поскольку волновая функция многомерная. Лучше всё-таки что сделано освоить. Допустим, симметрия сокращает вычисления. Но я в этом не спец.

По ответам автора темы видно, что тут не кванты, а ньютоновская механика.

Есть подозрение, что под "центрально-разностной" имеется в виду схема с перешагивание (leap-frog):

$v(t+\tau)=v(t-\tau)+\drac{F(t)}{m}\cdot\tau,$
$x(t+2\tau)=x(t)+v(t+\tau)\cdot\tau.$

[MathKvant]

Pphantom в сообщении #972222 писал(а):

На форуме принято обращение на "Вы". И чем Вам односторонняя аппроксимация не нравится?

Простите, не хотел обижать вас. Я попробовал реализовать одностороннюю.

Используется синтаксис Fortran

    
program molecul
    implicit none
    real*16, parameter :: r1=5.2917720859e-11! радиус первой орбиты в метрах
    real*16, parameter :: r2=21.1670883436e-11! радиус второй орбиты в метрах
    real*16, parameter :: me=9.10938291e-31! масса электронаб в кг
    real*16, parameter :: q=1.60217657e-19! заярд елетрона в кулонах
    real*16, parameter :: tau=1.52e-18! шаг по времени в секундах
    real*16 v_next_x1,v_previous_x1,x_next1,x_previous1,v_next_y1,v_previous_y1,y_next1,y_previous1
    real*16 v_next_x2,v_previous_x2,x_next2,x_previous2,v_next_y2,v_previous_y2,y_next2,y_previous2
    real*16 v_next_x3,v_previous_x3,x_next3,x_previous3,v_next_y3,v_previous_y3,y_next3,y_previous3
    real*16 x11,x22,x33,y11,y22,y33
    integer i
    x11=r1! начальные условия
    x22=r1
    x33=r2
    y11=0
    y22=0
    y33=0
        v_previous_x1=0
        x_previous1=0
        v_previous_y1=0
        y_previous1=0


        v_previous_x2=0
        x_previous2=0
        v_previous_y2=0
        y_previous2=0


        v_previous_x3=0
        x_previous3=0
        v_previous_y3=0
        y_previous3=0
    open(12,FILE='result.txt')
        do i=1,10000
        v_next_x1=v_previous_x1+tau*(q*q*abs(x22-x11)/((abs(x22-x11))**2+(abs(y22-y11))**2)+q*q*abs(x33-x11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-q*q*abs(-x11)/((abs(-x11))**2+(abs(-y11))**2))/m
        x_next1=x_previous1+tau*v_previous_x1
        v_next_y1=v_previous_y1+tau*(q*q*abs(y22-y11)/((abs(x22-x11))**2+(abs(y22-y11))**2)+q*q*abs(y33-y11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-q*q*abs(-y11)/((abs(-x11))**2+(abs(-y11))**2))/m
                y_next1=y_previous1+tau*v_previous_y1

        v_next_x2=v_previous_x2+tau*(q*q*abs(x11-x22)/((abs(x11-x22))**2+(abs(y11-y22))**2)+q*q*abs(x33-x22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-q*q*abs(-x22)/((abs(-x22))**2+(abs(-y22))**2))/m
        x_next2=x_previous2+tau*v_previous_x2
        v_next_y2=v_previous_y2+tau*(q*q*abs(y11-y22)/((abs(x11-x22))**2+(abs(y11-y22))**2)+q*q*abs(y33-y22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-q*q*abs(-y22)/((abs(-x22))**2+(abs(-y22))**2))/m
                y_next2=y_previous2+tau*v_previous_y2

        v_next_x3=v_previous_x+tau*(q*q*abs(x11-x33)/((abs(x11-x33))**2+(abs(y11-y33))**2)+q*q*abs(x22-x33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-q*q*abs(-x33)/((abs(-x33))**2+(abs(-y33))**2))/m
        x_next3=x_previous3+tau*v_previous_x3
        v_next_y3=v_previous_y3+tau*(q*q*abs(y11-y33)/((abs(x11-x33))**2+(abs(y11-y33))**2)+q*q*abs(y22-y33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-q*q*abs(-y33)/((abs(-x33))**2+(abs(-y33))**2))/m
                y_next3=y_previous3+tau*v_previous_y3

        write(12,*) x_next1,y_next1,x_next2,y_next2,x_next3,y_next3
        x11=x_next1
        x22=x_next2
        x33=x_next3
        y11=y_next1
        y22=y_next2
        y33=y_next3

        v_previous_x1=v_next_x1
        x_previous1=x_next1
        v_previous_y1=v_next_y1
        y_previous1=y_next1


        v_previous_x2=v_next_x2
        x_previous2=x_next
        v_previous_y2=v_next_y2
        y_previous2=y_next2


        v_previous_x3=v_next_x3
        x_previous3=x_next3
        v_previous_y3=v_next_y3
        y_previous3=y_next3
        enddo
    close(12)
    end program
 

Мне главное понять, верно ли я понимаю алгоритм. Программа, если что не работает, так как почему-то ругается на функцию abs().

DimaM в сообщении #972449 писал(а):

Есть подозрение, что под "центрально-разностной" имеется в виду схема с перешагивание (leap-frog):

Спасибо, похоже Вы правы. Теперь можно нормально гуглить и смотреть что к чему

[Pphantom]

DimaM в сообщении #972449 писал(а):

Есть подозрение, что под "центрально-разностной" имеется в виду схема с перешагивание (leap-frog):

Может быть, но ТС в этом не признается.

DimaM в сообщении #972449 писал(а):

$v(t+\tau)=v(t-\tau)+\drac{F(t)}{m}\cdot\tau,$
$x(t+2\tau)=x(t)+v(t+\tau)\cdot\tau.$

У Вас тут, похоже, три ошибки.

MathKvant в сообщении #972464 писал(а):

Я попробовал реализовать одностороннюю.

Как говорила Мальвина, "какой ужас!"

Конкретная проблема с функцией abs сводится к тому, что компиляторы Фортрана по умолчанию ограничивают длину строки кода 132 символами, а строчки, в которых Вы получаете ошибку, длиннее. abs просто не повезло - обрезка происходит как раз там, где начинается аргумент. Соответственно, надо либо переносить строчку (в какой-то момент ставите & и продолжаете на следующей), либо ключом компилятора увеличивать допустимую длину строки.

Но в целом код действительно ужасен. Во-первых, зачем Вам real*16? Для всего этого двойной точности более чем достаточно. Во-вторых, считать все в СИ, без предварительного обезразмеривания задачи, не стоит - это мощный источник как машинных погрешностей, так и ошибок при написании кода. Наконец, использовать Фортран и не пользоваться при этом векторными операциями - это реализация анекдота "взял билет и пошел пешком". Ну а шедевры вроде abs(-x) - это уже просто нет слов.

MathKvant в сообщении #972464 писал(а):

Мне главное понять, верно ли я понимаю алгоритм.

Боюсь, что в данном случае главная Ваша проблема связана с технологией программирования, а не с пониманием алгоритма. Если Вы займетесь реализацией leap-frog в том же стиле, результат, скорее всего, будет еще более жутким.

[DimaM]

Pphantom в сообщении #972502 писал(а):

У Вас тут, похоже, три ошибки.

Действительно :(.

[MathKvant]

Pphantom в сообщении #972502 писал(а):

Во-первых, зачем Вам real*16?

Да, действительно real*8 хватит, ставил с запасом.

Pphantom в сообщении #972502 писал(а):

Во-вторых, считать все в СИ, без предварительного обезразмеривания задачи, не стоит - это мощный источник как машинных погрешностей, так и ошибок при написании кода.

Здесь тоже вы правы, не могли подсказать, как лучше обезразмерить, в какие единицы перевести?

Pphantom в сообщении #972502 писал(а):

Наконец, использовать Фортран и не пользоваться при этом векторными операциями - это реализация анекдота "взял билет и пошел пешком"

Если честно, я не знаю, что за векторные операции и как ими пользоваться, если вы посоветуете, что почитать, буду очень благодарен.

Pphantom в сообщении #972502 писал(а):

Ну а шедевры вроде abs(-x) - это уже просто нет слов.

Я писал все в координатной форме, чтобы не запутаться, поэтому некоторые моменты нелепы.
Спасибо большое за конструктивную критику.

[Pphantom]

MathKvant в сообщении #973108 писал(а):

Здесь тоже вы правы, не могли подсказать, как лучше обезразмерить, в какие единицы перевести?

Если записать уравнения движения каждого электрона (движением ядра, насколько можно судить по коду, Вы пренебрегаете, что с одной стороны логично, с другой - на самом деле позволяет ограничиться даже одинарной точностью, т.е. real*4, при счете), окажется, что во все выражения входит одна и та же комбинация размерных постоянных $k\,e^2/m$, где $e$ - элементарный заряд, $m$ - масса электрона, $k$ - коэффициент в законе Кулона (о котором, кстати, Вы начисто забыли). Целесообразно выбрать систему единиц, в которой этот коэффициент единичен. Если при этом в качестве пространственного масштаба выбрать характерные величины для атома (например, $1~\mathring{A} = 10^{-10}$ м, то при этом у Вас получится и удобная единица времени.

MathKvant в сообщении #973108 писал(а):

Если честно, я не знаю, что за векторные операции и как ими пользоваться, если вы посоветуете, что почитать, буду очень благодарен.

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

MathKvant в сообщении #973108 писал(а):

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

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

[Munin]

Ландау, Лифшиц. Квантовая механика.



<...>

(Оффтоп)

Примечание: 2 Ry имеет тоже имя собственное: хартри, Ha. Хотя оно менее известно.

[MathKvant]

Pphantom в сообщении #973122 писал(а):

Если записать уравнения движения каждого электрона (движением ядра, насколько можно судить по коду, Вы пренебрегаете, что с одной стороны логично, с другой - на самом деле позволяет ограничиться даже одинарной точностью, т.е. real*4, при счете), окажется, что во все выражения входит одна и та же комбинация размерных постоянных $k\,e^2/m$, где $e$ - элементарный заряд, $m$ - масса электрона, $k$ - коэффициент в законе Кулона (о котором, кстати, Вы начисто забыли). Целесообразно выбрать систему единиц, в которой этот коэффициент единичен. Если при этом в качестве пространственного масштаба выбрать характерные величины для атома (например, $1~\mathring{A} = 10^{-10}$ м, то при этом у Вас получится и удобная единица времени.

Спасибо, я вроде подправил самые тупые моменты,перевел все в СГСЭ,теперь коэф k=1,также вынес степень $10^{-10}$, теперь все компиллится, но значения все NaN.

Используется синтаксис Fortran

    program molecul
    implicit none
    real*8, parameter :: r1=0.52917720859! радиус первой орбиты в ангстрем
    real*8, parameter :: r2=2.11670883436! радиус второй орбиты в ангстрем
    real*8, parameter :: me=9.10938291e-18! масса электрона в гр
    real*8, parameter :: q=4.8032042713! заряд элетрона в сгсэ
    real*8, parameter :: tau=1.52e-8! шаг по времени в секундах
    real*8 v_next_x1,v_previous_x1,x_next1,x_previous1,v_next_y1,v_previous_y1,y_next1,y_previous1
    real*8 v_next_x2,v_previous_x2,x_next2,x_previous2,v_next_y2,v_previous_y2,y_next2,y_previous2
    real*8 v_next_x3,v_previous_x3,x_next3,x_previous3,v_next_y3,v_previous_y3,y_next3,y_previous3
    real*8 x11,x22,x33,y11,y22,y33
    real*8 const
    integer i
    const=tau*q*q/me
    write(*,*) const
    x11=r1! начальные условия
    x22=r1
    x33=r2
    y11=0
    y22=0
    y33=0
        v_previous_x1=0
        x_previous1=0
        v_previous_y1=0
        y_previous1=0


        v_previous_x2=0
        x_previous2=0
        v_previous_y2=0
        y_previous2=0


        v_previous_x3=0
        x_previous3=0
        v_previous_y3=0
        y_previous3=0
    open(12,FILE='result.txt')
        do i=1,5
        v_next_x1=v_previous_x1+tau*(q*q*abs(x22-x11)/((abs(x22-x11))**2+(abs(y22-y11))**2)&
        +q*q*abs(x33-x11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-q*q*abs(x11)/((abs(x11))**2+(abs(y11))**2))/me
        x_next1=x_previous1+tau*v_previous_x1
        v_next_y1=v_previous_y1+tau*(q*q*abs(y22-y11)/((abs(x22-x11))**2+(abs(y22-y11))**2)&
        +q*q*abs(y33-y11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-q*q*abs(y11)/((abs(x11))**2+(abs(y11))**2))/me
                y_next1=y_previous1+tau*v_previous_y1

        v_next_x2=v_previous_x2+tau*(q*q*abs(x11-x22)/((abs(x11-x22))**2+(abs(y11-y22))**2)&
        +q*q*abs(x33-x22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-q*q*abs(x22)/((abs(x22))**2+(abs(y22))**2))/me
        x_next2=x_previous2+tau*v_previous_x2
        v_next_y2=v_previous_y2+tau*(q*q*abs(y11-y22)/((abs(x11-x22))**2+(abs(y11-y22))**2)&
        +q*q*abs(y33-y22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-q*q*abs(y22)/((abs(x22))**2+(abs(y22))**2))/me
                y_next2=y_previous2+tau*v_previous_y2

        v_next_x3=v_previous_x3+tau*(q*q*abs(x11-x33)/((abs(x11-x33))**2+(abs(y11-y33))**2)&
        +q*q*abs(x22-x33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-q*q*abs(x33)/((abs(x33))**2+(abs(y33))**2))/me
        x_next3=x_previous3+tau*v_previous_x3
        v_next_y3=v_previous_y3+tau*(q*q*abs(y11-y33)/((abs(x11-x33))**2+(abs(y11-y33))**2)&
        +q*q*abs(y22-y33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-q*q*abs(y33)/((abs(x33))**2+(abs(y33))**2))/me
                y_next3=y_previous3+tau*v_previous_y3

        write(12,*) x_next1,y_next1,x_next2,y_next2,x_next3,y_next3
        x11=x_next1
        x22=x_next2
        x33=x_next3
        y11=y_next1
        y22=y_next2
        y33=y_next3

        v_previous_x1=v_next_x1
        x_previous1=x_next1
        v_previous_y1=v_next_y1
        y_previous1=y_next1


        v_previous_x2=v_next_x2
        x_previous2=x_next2
        v_previous_y2=v_next_y2
        y_previous2=y_next2


        v_previous_x3=v_next_x3
        x_previous3=x_next3
        v_previous_y3=v_next_y3
        y_previous3=y_next3
        enddo
    close(12)
    end program
 

Munin, спасибо!