2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 14:25 


22/09/10
75
Привет всем, возник вопрос по МД. Задача состоит в нахождении траектории движения для каждого электрона в атоме лития. Для нахождения скорости и положения в пространстве использую центрально-разностную апроксимацию.
В итоге имею вот такие уравнения:
$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 нужно подставлять для начала подсчета.

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 14:45 
Заслуженный участник


09/05/12
25179
MathKvant в сообщении #972193 писал(а):
Привет всем, возник вопрос по МД. Задача состоит в нахождении траектории движения для каждого электрона в атоме лития.
Во-первых, что такое "МД"? Молекулярная динамика?

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

MathKvant в сообщении #972193 писал(а):
Для нахождения скорости и положения в пространстве использую центрально-разностную апроксимацию.
И как, получилось? :wink: Собственно, первый вопрос, который тут появляется - зачем Вы ее используете?

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 14:59 


22/09/10
75
Pphantom в сообщении #972201 писал(а):
Во-первых, что такое "МД"? Молекулярная динамика?

Да, оно самое
Pphantom в сообщении #972201 писал(а):
Во-вторых, задача в такой постановке, мягко говоря, выглядит несколько странно. Вы уверены, что у электрона в атоме есть траектория?

Не я придумал задачу, может и глупо рассматривать траектория движения электрона, я это представляю как очень маленькие планеты.
Pphantom в сообщении #972201 писал(а):
И как, получилось? :wink: Собственно, первый вопрос, который тут появляется - зачем Вы ее используете?

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 15:08 
Заслуженный участник


09/05/12
25179
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}
$
Вопрос не в том, что такое центрально-разностная аппроксимация (я это знаю), а в том, зачем Вы ее используете. Просто если Вы пытаетесь таким образом осложнить себе жизнь во имя какой-то великой цели :D , то это одно, а если это было сделано "просто так", то проще просто так не делать.

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 15:15 


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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 15:37 
Заслуженный участник


09/05/12
25179
MathKvant в сообщении #972215 писал(а):
Подскажи ,пожалуйста, как тогда правильно решать мою задачу.
На форуме принято обращение на "Вы". И чем Вам односторонняя аппроксимация не нравится?

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение01.02.2015, 17:13 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение02.02.2015, 10:33 
Заслуженный участник


28/12/12
7740
мат-ламер в сообщении #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.$

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение02.02.2015, 12:16 


22/09/10
75
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):

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение02.02.2015, 14:13 
Заслуженный участник


09/05/12
25179
DimaM в сообщении #972449 писал(а):
Есть подозрение, что под "центрально-разностной" имеется в виду схема с перешагивание (leap-frog):
Может быть, но ТС в этом не признается. :D

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 писал(а):
Я попробовал реализовать одностороннюю.
Как говорила Мальвина, "какой ужас!" :facepalm:

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

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

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение02.02.2015, 17:09 
Заслуженный участник


28/12/12
7740
Pphantom в сообщении #972502 писал(а):
У Вас тут, похоже, три ошибки.

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение03.02.2015, 18:05 


22/09/10
75
Pphantom в сообщении #972502 писал(а):
Во-первых, зачем Вам real*16?

Да, действительно real*8 хватит, ставил с запасом.
Pphantom в сообщении #972502 писал(а):
Во-вторых, считать все в СИ, без предварительного обезразмеривания задачи, не стоит - это мощный источник как машинных погрешностей, так и ошибок при написании кода.

Здесь тоже вы правы, не могли подсказать, как лучше обезразмерить, в какие единицы перевести?
Pphantom в сообщении #972502 писал(а):
Наконец, использовать Фортран и не пользоваться при этом векторными операциями - это реализация анекдота "взял билет и пошел пешком"

Если честно, я не знаю, что за векторные операции и как ими пользоваться, если вы посоветуете, что почитать, буду очень благодарен.
Pphantom в сообщении #972502 писал(а):
Ну а шедевры вроде abs(-x) - это уже просто нет слов.

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение03.02.2015, 18:46 
Заслуженный участник


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

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

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение03.02.2015, 19:37 
Заслуженный участник
Аватара пользователя


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

Изображение

<...>

Изображение
Изображение

(Оффтоп)

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

 Профиль  
                  
 
 Re: расчет движения электронов в атоме лития методом МД(2мерная)
Сообщение06.02.2015, 12:11 


22/09/10
75
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, спасибо!

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 24 ]  На страницу 1, 2  След.

Модераторы: photon, whiterussian, profrotter, Jnrty, Aer, Парджеттер, Eule_A, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: lazarius


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group