Форум программистов, компьютерный форум CyberForum.ru

Метод Монте-Карло для систем многих частиц - C++

Войти
Регистрация
Восстановить пароль
 
hellmoore
0 / 0 / 0
Регистрация: 07.04.2015
Сообщений: 1
11.04.2016, 13:06     Метод Монте-Карло для систем многих частиц #1
Есть задача реализовать метод Монте-Карло для систем многих частиц.

На данный момент сделан расчет взаимодействия двух частиц, а также граничные условия.

Далее надо 1) реализовать схему начальной расстановки частиц(кристалл или случайное распределение), проверить корректность работы для различного числа частиц; 2) Реализовать схему Метрополиса. Подобрать амплитуду случайной силы.

Проблема в том, что нет особого понимания, так как с физикой не сталкивалась со школы, а как это запрограммировать вообще не знаю.

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <iomanip>
 
using namespace std;
 
const int ntot = 20000;
FILE *out = stdout;
 
 
class MD
{
    private:
        const int N = 5; // количество частиц
        double dt;
    double rn[N][3], L[3], L2[3];
        double r[N][3];
        double f[N][3]; // координаты, сила
        double utot; 
        double E_kin;
    public:
        double v[N][3];
 
        MD()
        {
            N = 5;
            dt = 0.0005;
            utot = 0;
            E_kin = 0;
 
            for(int i = 0; i < 3; i++)
            {
                L[i] = 3;
                L2[i] = L[i] / 2;
            }
 
            out = fopen("result.csv","ab+");
        }
 
 
    void nearest_image(); // приведение к основной ячейке
        void ClearForces();// f = 0;
        void ReadCoordinates(); // чтение координат из файла
        void CalcForces(); //Расчет сил
        void EqMotion();
        void KineticEnergy();
        void SaveState(int n);
};
 
void MD::ClearForces()
{
 
    for (int i = 0; i < N; i ++)
        for (int j = 0; j < 3; j ++)
            f[i][j] =0;
}
 
void MD::ReadCoordinates()
{
    FILE *a = stdout;
    a = fopen("coordinates.csv","r");
 
    for (int i = 0; i < N; i ++)
        for (int j = 0; j < 3; j ++)
            fscanf(a,"%lf;", &r[i][j]);
 
    fclose(a);
}
 
void MD::CalcForces() //Расчет сил
{
 
 
    double rij[3]; // расстояние между частицами
    double Force = 0; // сила
    double r1, u1, u6, u12;
    utot = 0;
 
    for(int i = 1; i < 5; i ++)
        for(int j = 0; j < i; j ++)
        {
            double r2 = 0;
            for(int k=0; k<3; k++)
            {
                rij[k] = rn[i][k] - rn[j][k];
                if(rij[k] > L2[k]) rij[k] -= L[k];
                else if(rij[k] < -L2[k]) rij[k] += L[k];
                r2 += rij[k] * rij[k];
            }
 
            r1 = sqrt(r2);
            u1 = r1*r1*r1;
            u6 = u1*u1;
            u12 = u6*u6;
            utot += 4.*(1./u12 - 1./u6); // расчет потенциала
            Force = -24.*(u6-2.)/(u12 *r1*r1);
 
            for(int k = 0; k < 3; k ++)
            {
                f[i][k] += Force * rij[k];
                f[j][k] -= Force * rij[k];
            }
        }
}
 
void MD::EqMotion() //Решение ур.движения
{
    for(int i = 0; i < 5; i ++)
        for(int k = 0; k < 3; k ++)
        {
            v[i][k] += f[i][k] * dt;
            r[i][k] += v[i][k] * dt;
        }
}
 
void MD::KineticEnergy()
{
    double velocity = 0;
 
    for(int i = 0; i < 5; i ++)
        for(int j = 0; j < 3; j ++)
            velocity += v[i][j]*v[i][j];
 
        E_kin = velocity/2.;
}
 
void MD::nearest_image()
{
    for(int i=0; i<5; i++)
        for(int k=0; k<3; k++)
        {
            if(r[i][k] > 0)
                rn[i][k] = fmod(r[i][k] + L2[k], L[k]) - L2[k];
            else
                rn[i][k] = fmod(r[i][k] - L2[k], L[k]) + L2[k];
        }
 
}
 
void MD::SaveState(int n) // сохранение результата
{
    fprintf(out, "%.3lf; %.10lf; %.10lf; %.10lf;\n", n*dt, utot, E_kin, utot+E_kin);
}
 
int main()
{
 
    MD A;
    for (int i = 0; i < 5; i ++)
        for (int j = 0; j < 3; j ++)
            A.v[i][j] = 0;
 
    A.ReadCoordinates();
    
 
    for(int n = 0; n < ntot; n ++)
    {
    A.nearest_image();
        A.ClearForces();
        A.CalcForces();
        A.KineticEnergy();
        A.EqMotion();
 
        if(n%10 == 0) A.SaveState(n);
    }
    fclose(out);
    return 0;
}
Файл с координатами https://yadi.sk/i/8m7Sl2e3qodH3
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
11.04.2016, 13:06     Метод Монте-Карло для систем многих частиц
Посмотрите здесь:

метод Монте-Карло C++
C++ метод Монте-Карло
C++ Метод Монте-Карло
C++ Метод Монте-Карло(непонятная неработоспособность программы)
C++ Метод монте Карло
C++ вроде метод монте карло
Метод Монте-Карло. Объем сферы C++
C++ Моделирование систем массового обслуживания метод монте карло
C++ Метод Монте-Карло
C++ Метод Монте Карло (неправильные значения)
Метод Монте-Карло C++
C++ Метод Монте-Карло в вычислении площади многоугольника

Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
Ответ Создать тему

Метки
физика
Опции темы

Текущее время: 01:17. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Рейтинг@Mail.ru