-
Notifications
You must be signed in to change notification settings - Fork 0
/
ST_mod.f90
222 lines (187 loc) · 10.3 KB
/
ST_mod.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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
module ST_mod !Author: Miguel Benito de Lama (January 2023)
use general_mod
implicit none
contains
subroutine ST_potential_diff(ini_coor, trial_coor, hL, L, r_cut, N, atom, delta_v)
!Calculates the Stillinger potential energy difference of two given structures
implicit none
integer(kind=precision),intent(in):: N, atom
integer(kind=precision):: i, j, k
real(kind=precision),intent(in)::ini_coor(:,:), trial_coor(:,:), hL, L, r_cut
real(kind=precision):: result, r_atom_i_2,r_atom_k_2, r_i_k_2, theta_jik, theta_ijk
real(kind=precision),intent(out)::delta_v
delta_v = 0
do i=1,N
if (i .ne. atom) then
call distance_squared(ini_coor(:,atom),ini_coor(:,i), r_atom_i_2, hL, L)
if (r_atom_i_2 .lt. r_cut**2.) then
call f_2(r_atom_i_2,r_cut,result)
delta_v=delta_v-result
do k=i+1,N
if (k .ne. atom) then
call distance_squared(ini_coor(:,atom),ini_coor(:,k), r_atom_k_2, hL, L)
if (r_atom_k_2 .lt. r_cut**2.) then
call angle_jik(ini_coor(:,atom),ini_coor(:,i), ini_coor(:,k), hL, L, theta_jik) ! angle theta_jik
call h_f_3(r_atom_i_2, r_atom_k_2, r_cut, theta_jik, result)
delta_v=delta_v-result
endif
endif
enddo
do k=1,N
if ((k .ne. atom) .and. (k .ne. i)) then
call distance_squared(ini_coor(:,i),ini_coor(:,k), r_i_k_2, hL, L) ! distance rjk
if (r_i_k_2 .lt. r_cut**2.) then
call angle_jik(ini_coor(:,i),ini_coor(:,atom), ini_coor(:,k), hL, L, theta_ijk) ! angle theta_ijk
call h_f_3(r_atom_i_2, r_i_k_2, r_cut, theta_ijk, result)
delta_v=delta_v-result
endif
endif
enddo
endif
call distance_squared(trial_coor(:,atom),trial_coor(:,i), r_atom_i_2, hL, L)
if (r_atom_i_2 .lt. r_cut**2.) then
call f_2(r_atom_i_2,r_cut,result)
delta_v=delta_v+result
do k=i+1,N
if (k .ne. atom) then
call distance_squared(trial_coor(:,atom),trial_coor(:,k), r_atom_k_2, hL, L)
if (r_atom_k_2 .lt. r_cut**2.) then
call angle_jik(trial_coor(:,atom),trial_coor(:,i), trial_coor(:,k), hL, L, theta_jik) ! angle theta_jik
call h_f_3(r_atom_i_2, r_atom_k_2, r_cut, theta_jik, result)
delta_v=delta_v+result
endif
endif
enddo
do k=1,N
if ((k .ne. atom) .and. (k .ne. i)) then
call distance_squared(trial_coor(:,i),trial_coor(:,k), r_i_k_2, hL, L) ! distance rjk
if (r_i_k_2 .lt. r_cut**2.) then
call angle_jik(trial_coor(:,i),trial_coor(:,atom), trial_coor(:,k), hL, L, theta_ijk) ! angle theta_ijk
call h_f_3(r_atom_i_2, r_i_k_2, r_cut, theta_ijk, result)
delta_v=delta_v+result
endif
endif
enddo
endif
endif
enddo
end subroutine ST_potential_diff
subroutine f_2(r_ij_2,r_cut,result)
!Calculates the f_2 function used to obtain the V2 term of the potential
implicit none
real(kind=precision), intent(in)::r_cut, r_ij_2
real(kind=precision), parameter:: A=7.049556277, B=0.6022245584
real(kind=precision),intent(out):: result
result=A* ((B * r_ij_2**(-2.)) - 1) * EXP(1 / (sqrt(r_ij_2) - r_cut))
end subroutine f_2
subroutine angle_jik(atom_i, atom_j, atom_k, hL,L, theta)
!Calculates the angle jik between three given atoms, atom_i, atom_j and atom_k
implicit none
real(kind=precision), intent(in):: atom_i(3,1), atom_j(3,1), atom_k(3,1),hL,L
real(kind=precision)::prod_v_ij_v_ik, mod_v_ij, mod_v_ik, vec_ij(3), vec_ik(3)
real(kind=precision),intent(out):: theta
integer(kind=precision)::i
prod_v_ij_v_ik=0
mod_v_ij=0
mod_v_ik=0
do i=1,3
vec_ij(i)=atom_j(i,1)-atom_i(i,1)
if (vec_ij(i) .gt. hL) then
vec_ij(i)=vec_ij(i)-L
elseif(vec_ij(i) .lt. -hL) then
vec_ij(i)=vec_ij(i)+L
endif
vec_ik(i)=atom_k(i,1)-atom_i(i,1)
if (vec_ik(i) .gt. hL) then
vec_ik(i)=vec_ik(i)-L
elseif(vec_ik(i) .lt. -hL) then
vec_ik(i)=vec_ik(i)+L
endif
mod_v_ij=mod_v_ij+vec_ij(i)**(2.)
mod_v_ik=mod_v_ik+vec_ik(i)**(2.)
prod_v_ij_v_ik=prod_v_ij_v_ik+vec_ij(i)*vec_ik(i)
enddo
theta=ACOS((prod_v_ij_v_ik)/(sqrt(mod_v_ij)*sqrt(mod_v_ik)))
end subroutine angle_jik
subroutine h_f_3(r_ij_2, r_ik_2, r_cut, theta, result)
!Calculates the h function used to obtain the V3 term of the potential
implicit none
real(kind=precision), intent(in):: r_ij_2, r_ik_2, r_cut, theta
real(kind=precision), parameter:: lambda=21.0, gamma=1.20
real(kind=precision):: exp_term, cos_term
real(kind=precision),intent(out):: result
exp_term = (gamma/(sqrt(r_ij_2)-r_cut))+(gamma/(sqrt(r_ik_2)-r_cut))
cos_term = (COS(theta)+(1./3.))**(2.)
result=lambda *EXP(exp_term) * cos_term
end subroutine h_f_3
subroutine ST_potential_diff_NL(ini_coor, trial_coor, neighbour_matrix, hL, L, r_cut, atom, n_max, delta_v)
!Calculates the Stillinger potential energy difference of two given structures using neighbour lists
implicit none
integer(kind=precision),intent(in):: atom,n_max
integer,intent(in)::neighbour_matrix(:,:)
integer(kind=precision):: i, j, k
real(kind=precision),intent(in)::ini_coor(:,:), trial_coor(:,:), hL, L, r_cut
real(kind=precision)::f_2_result, h_f_3_result,r_atom_i_2,r_atom_k_2, r_i_k_2, theta_jik, theta_ijk
real(kind=precision),intent(out)::delta_v
delta_v = 0
do i=1,n_max
if (neighbour_matrix(i,atom).eq. 0) exit
call distance_squared(ini_coor(:,atom),ini_coor(:,neighbour_matrix(i,atom)), r_atom_i_2, hL, L)
if (r_atom_i_2 .lt. r_cut**2.) then
call f_2(r_atom_i_2,r_cut,f_2_result)
delta_v=delta_v-f_2_result !V2 term
do k=i+1,n_max
if (neighbour_matrix(k,atom).eq. 0) exit
call distance_squared(ini_coor(:,atom),ini_coor(:,neighbour_matrix(k,atom)), r_atom_k_2, hL, L)
if (r_atom_k_2 .lt. r_cut**2.) then
call angle_jik(ini_coor(:,atom),ini_coor(:,neighbour_matrix(i,atom)), &
ini_coor(:,neighbour_matrix(k,atom)), hL, L, theta_jik)
call h_f_3(r_atom_i_2, r_atom_k_2, r_cut, theta_jik, h_f_3_result)
delta_v=delta_v-h_f_3_result !V3 term
endif
enddo
do k=1,n_max
if (neighbour_matrix(k,neighbour_matrix(i,atom)).eq.0) exit
if (neighbour_matrix(k,neighbour_matrix(i,atom)) .ne. atom) then
call distance_squared(ini_coor(:,neighbour_matrix(i,atom)),&
ini_coor(:,neighbour_matrix(k,neighbour_matrix(i,atom))), r_i_k_2, hL, L)
if (r_i_k_2 .lt. r_cut**2.) then
call angle_jik(ini_coor(:,neighbour_matrix(i,atom)),ini_coor(:,atom), &
ini_coor(:,neighbour_matrix(k,neighbour_matrix(i,atom))), hL, L, theta_ijk)
call h_f_3(r_atom_i_2, r_i_k_2, r_cut, theta_ijk, h_f_3_result)
delta_v=delta_v-h_f_3_result !V3 term
endif
endif
enddo
endif
call distance_squared(trial_coor(:,atom),trial_coor(:,neighbour_matrix(i,atom)), r_atom_i_2, hL, L)
if (r_atom_i_2 .lt. r_cut**2.) then
call f_2(r_atom_i_2,r_cut,f_2_result)
delta_v=delta_v+f_2_result !V2 term
do k=i+1,n_max
if (neighbour_matrix(k,atom).eq. 0) exit
call distance_squared(trial_coor(:,atom),trial_coor(:,neighbour_matrix(k,atom)), r_atom_k_2, hL, L)
if (r_atom_k_2 .lt. r_cut**2.) then
call angle_jik(trial_coor(:,atom),trial_coor(:,neighbour_matrix(i,atom)),&
trial_coor(:,neighbour_matrix(k,atom)), hL, L, theta_jik)
call h_f_3(r_atom_i_2, r_atom_k_2, r_cut, theta_jik, h_f_3_result)
delta_v=delta_v+h_f_3_result !V3 term
endif
enddo
do k=1,n_max
if (neighbour_matrix(k,neighbour_matrix(i,atom)).eq.0) exit
if (neighbour_matrix(k,neighbour_matrix(i,atom)) .ne. atom) then
call distance_squared(trial_coor(:,neighbour_matrix(i,atom)),&
trial_coor(:,neighbour_matrix(k,neighbour_matrix(i,atom))), r_i_k_2, hL, L)
if (r_i_k_2 .lt. r_cut**2.) then
call angle_jik(trial_coor(:,neighbour_matrix(i,atom)),trial_coor(:,atom), &
trial_coor(:,neighbour_matrix(k,neighbour_matrix(i,atom))), hL, L, theta_ijk)
call h_f_3(r_atom_i_2, r_i_k_2, r_cut, theta_ijk, h_f_3_result)
delta_v=delta_v+h_f_3_result !V3 term
endif
endif
enddo
endif
enddo
end subroutine ST_potential_diff_NL
end module ST_mod