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
6594
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, Супермодераторы



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

Сейчас этот форум просматривают: нет зарегистрированных пользователей


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

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