-
Notifications
You must be signed in to change notification settings - Fork 20
/
Copy pathexample_LJ_min_image_model.F90
164 lines (136 loc) · 4.45 KB
/
example_LJ_min_image_model.F90
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
! API documented in example_LJ_model.F90
subroutine ll_init_model(N_params, params)
use example_LJ_params_mod
implicit none
integer :: N_params
double precision :: params(N_params)
epsilon(1,1) = 4.0
epsilon(2,2) = 4.0
epsilon(1,2) = 6.0
epsilon(2,1) = epsilon(1,2)
sigma(1,1) = 1.0
sigma(2,2) = 1.0
sigma(1,2) = (sigma(1,1)+sigma(2,2))/2.0
sigma(2,1) = sigma(1,2)
cutoff = 3.0*sigma
cutoff_sq = cutoff*cutoff
E_offset = (sigma/cutoff)**12 - (sigma/cutoff)**6
end subroutine ll_init_model
subroutine ll_init_config(N, Z, pos, cell, Emax)
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), cell(3,3)
double precision :: Emax
return
end subroutine ll_init_config
double precision function ll_eval_energy(N, Z, pos, n_extra_data, extra_data, cell)
use example_mat_mod
use example_LJ_params_mod
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), cell(3,3)
integer :: n_extra_data
double precision :: extra_data(n_extra_data, N)
integer :: i, j, Z_i, Z_j
double precision :: dr(3), dr_mag, dr_l(3)
double precision :: cell_inv(3,3)
call matrix3x3_inverse(cell, cell_inv)
ll_eval_energy = 0.0
do i=1, N
Z_i = Z(i)
do j=i+1, N
Z_j = Z(j)
dr = pos(:,i)-pos(:,j)
dr_l = matmul(cell_inv, dr)
dr_l = dr_l - floor(dr_l+0.5)
dr = matmul(cell, dr_l)
dr_mag = sqrt(sum(dr*dr))
if (dr_mag < cutoff(Z_i,Z_j)) then
ll_eval_energy = ll_eval_energy + epsilon(Z_i,Z_j)*((sigma(Z_i,Z_j)/dr_mag)**12 - &
(sigma(Z_i,Z_j)/dr_mag)**6 - E_offset(Z_i,Z_j))
endif
end do
end do
end function ll_eval_energy
integer function ll_move_atom_1(N, Z, pos, n_extra_data, extra_data, cell, d_i, d_pos, dEmax, dE)
use example_mat_mod
use example_LJ_params_mod
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), cell(3,3)
integer :: n_extra_data
double precision :: extra_data(n_extra_data, N)
integer :: d_i
double precision :: d_pos(3)
double precision :: dEmax, dE
integer :: i, j, Z_i, Z_j
double precision :: dr0(3), dr1(3), dr0_l(3), dr1_l(3), dr0_mag, dr1_mag
double precision :: cell_inv(3,3)
call matrix3x3_inverse(cell, cell_inv)
dE = 0.0
i=d_i
Z_i = Z(i)
do j=1,N
if (j == i) cycle
Z_j = Z(j)
dr0 = pos(:,i) - pos(:,j)
dr1 = (pos(:,i)+d_pos(:)) - pos(:,j)
dr0_l = matmul(cell_inv, dr0)
dr0_l = dr0_l - floor(dr0_l+0.5)
dr0 = matmul(cell, dr0_l)
dr0_mag = sqrt(sum(dr0*dr0))
dr1_l = matmul(cell_inv, dr1)
dr1_l = dr1_l - floor(dr1_l+0.5)
dr1 = matmul(cell, dr1_l)
dr1_mag = sqrt(sum(dr1*dr1))
if (dr0_mag < cutoff(Z_i,Z_j)) then
dE = dE - epsilon(Z_i,Z_j)*((sigma(Z_i,Z_j)/dr0_mag)**12 - (sigma(Z_i,Z_j)/dr0_mag)**6 - E_offset(Z_i,Z_j))
endif
if (dr1_mag < cutoff(Z_i,Z_j)) then
dE = dE + epsilon(Z_i,Z_j)*((sigma(Z_i,Z_j)/dr1_mag)**12 - (sigma(Z_i,Z_j)/dr1_mag)**6 - E_offset(Z_i,Z_j))
endif
end do
if (dE < dEmax) then ! accept
pos(:,i) = pos(:,i) + d_pos(:)
ll_move_atom_1 = 1
else ! reject
dE = 0.0
ll_move_atom_1 = 0
endif
end function ll_move_atom_1
function ll_eval_forces(N, Z, pos, n_extra_data, extra_data, cell, forces) result(energy)
use example_mat_mod
use example_LJ_params_mod
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), cell(3,3), forces(3,N)
integer :: n_extra_data
double precision :: extra_data(n_extra_data, N)
double precision :: energy ! result
integer :: i, j, Z_i, Z_j
double precision :: dr(3), dr_mag, dr_l(3)
double precision :: cell_inv(3,3)
call matrix3x3_inverse(cell, cell_inv)
energy = 0.0
forces = 0.0
do i=1, N
Z_i = Z(i)
do j=i+1, N
Z_j = Z(j)
dr = pos(:,i)-pos(:,j)
dr_l = matmul(cell_inv, dr)
dr_l = dr_l - floor(dr_l+0.5)
dr = matmul(cell, dr_l)
dr_mag = sqrt(sum(dr*dr))
if (dr_mag < cutoff(Z_i,Z_j)) then
energy = energy + epsilon(Z_i,Z_j)*((sigma(Z_i,Z_j)/dr_mag)**12 - (sigma(Z_i,Z_j)/dr_mag)**6 - E_offset(Z_i,Z_j))
forces(:,i) = forces(:,i) - epsilon(Z_i,Z_j)*(-12.0*sigma(Z_i,Z_j)**12/dr_mag**13 + 6.0*sigma(Z_i,Z_j)**6/dr_mag**7)*(dr/dr_mag)
forces(:,j) = forces(:,j) + epsilon(Z_i,Z_j)*(-12.0*sigma(Z_i,Z_j)**12/dr_mag**13 + 6.0*sigma(Z_i,Z_j)**6/dr_mag**7)*(dr/dr_mag)
endif
end do
end do
end function ll_eval_forces