1 |
chuckv |
4 |
module model_module |
2 |
|
|
use definitions, ONLY: DP,ndim,machdep_lnblnk |
3 |
|
|
implicit none |
4 |
|
|
PRIVATE |
5 |
|
|
|
6 |
|
|
include 'headers/atom.h' |
7 |
|
|
integer,public, parameter :: natypes = 202 |
8 |
|
|
|
9 |
|
|
public :: is_atype |
10 |
|
|
public :: atype_mass |
11 |
|
|
public :: atype_atoi |
12 |
|
|
public :: atype_identifier |
13 |
|
|
public :: get_max_atype |
14 |
|
|
public :: model_na |
15 |
|
|
public :: model_atype |
16 |
|
|
public :: set_pos |
17 |
|
|
! public :: get_max_atomindex |
18 |
|
|
|
19 |
|
|
contains |
20 |
|
|
!logical |
21 |
|
|
function is_atype(atype) result(junk) |
22 |
|
|
integer :: atype |
23 |
|
|
logical :: junk |
24 |
|
|
|
25 |
|
|
select case(atype) |
26 |
|
|
case(H_atom) |
27 |
|
|
junk = .true. |
28 |
|
|
case(He_atom) |
29 |
|
|
junk = .true. |
30 |
|
|
case(C_atom) |
31 |
|
|
junk = .true. |
32 |
|
|
case(N_atom) |
33 |
|
|
junk = .true. |
34 |
|
|
case(O_atom) |
35 |
|
|
junk = .true. |
36 |
|
|
case(F_atom) |
37 |
|
|
junk = .true. |
38 |
|
|
case(Ne_atom) |
39 |
|
|
junk = .true. |
40 |
|
|
case(S_atom) |
41 |
|
|
junk = .true. |
42 |
|
|
case(Cl_atom) |
43 |
|
|
junk = .true. |
44 |
|
|
case(Ar_atom) |
45 |
|
|
junk = .true. |
46 |
|
|
case(Ni_atom) |
47 |
|
|
junk = .true. |
48 |
|
|
case(Cu_atom) |
49 |
|
|
junk = .true. |
50 |
|
|
case(Rh_atom) |
51 |
|
|
junk = .true. |
52 |
|
|
case(Pd_atom) |
53 |
|
|
junk = .true. |
54 |
|
|
case(Ag_atom) |
55 |
|
|
junk = .true. |
56 |
|
|
case(Ir_atom) |
57 |
|
|
junk = .true. |
58 |
|
|
case(Pt_atom) |
59 |
|
|
junk = .true. |
60 |
|
|
case(Au_atom) |
61 |
|
|
junk = .true. |
62 |
|
|
case(Pb_atom) |
63 |
|
|
junk = .true. |
64 |
|
|
case(Al_atom) |
65 |
|
|
junk = .true. |
66 |
|
|
case(Br_atom) |
67 |
|
|
junk = .true. |
68 |
|
|
case(Kr_atom) |
69 |
|
|
junk = .true. |
70 |
|
|
case(HV_atom) |
71 |
|
|
junk = .true. |
72 |
|
|
case(LT_atom) |
73 |
|
|
junk = .true. |
74 |
|
|
case(DZ_atom) |
75 |
|
|
junk = .true. |
76 |
|
|
case default |
77 |
|
|
junk = .false. |
78 |
|
|
end select |
79 |
|
|
|
80 |
|
|
end function is_atype |
81 |
|
|
|
82 |
|
|
!real double returns 0.0E0_DP if error. |
83 |
|
|
function atype_mass(atype) result(m) |
84 |
|
|
|
85 |
|
|
integer, intent(in) :: atype |
86 |
|
|
real( kind = DP ) :: m |
87 |
|
|
|
88 |
|
|
select case(atype) |
89 |
|
|
case(H_atom) |
90 |
|
|
m = 1.00794E0_DP |
91 |
|
|
case(He_atom) |
92 |
|
|
m = 4.003E0_DP |
93 |
|
|
case(C_atom) |
94 |
|
|
m = 12.011E0_DP |
95 |
|
|
case(N_atom) |
96 |
|
|
m = 14.01E0_DP |
97 |
|
|
case(O_atom) |
98 |
|
|
m = 15.9994E0_DP |
99 |
|
|
case(F_atom) |
100 |
|
|
m = 19.00E0_DP |
101 |
|
|
case(Ne_atom) |
102 |
|
|
m = 20.18E0_DP |
103 |
|
|
case(S_atom) |
104 |
|
|
m = 32.066E0_DP |
105 |
|
|
case(Cl_atom) |
106 |
|
|
m = 35.45E0_DP |
107 |
|
|
case(Ar_atom) |
108 |
|
|
m = 39.948E0_DP |
109 |
|
|
case(Ag_atom) |
110 |
|
|
m = 107.868E0_DP |
111 |
|
|
case(Cu_atom) |
112 |
|
|
m = 63.546E0_DP |
113 |
|
|
case(Ni_atom) |
114 |
|
|
m = 58.691E0_DP |
115 |
|
|
case(Rh_atom) |
116 |
|
|
m = 102.90550E0_DP |
117 |
|
|
case(Pd_atom) |
118 |
|
|
m = 106.42E0_DP |
119 |
|
|
case(Ir_atom) |
120 |
|
|
m = 192.22E0_DP |
121 |
|
|
case(Pt_atom) |
122 |
|
|
m = 195.08E0_DP |
123 |
|
|
case(Au_atom) |
124 |
|
|
m = 196.96654E0_DP |
125 |
|
|
case(Pb_atom) |
126 |
|
|
m = 207.2E0_DP |
127 |
|
|
case(Al_atom) |
128 |
|
|
m = 26.981539E0_DP |
129 |
|
|
case(Br_atom) |
130 |
|
|
m = 79.91E0_DP |
131 |
|
|
case(Kr_atom) |
132 |
|
|
m = 83.80E0_DP |
133 |
|
|
case(HV_atom) |
134 |
|
|
m = 100.0E0_DP * 39.948E0_DP |
135 |
|
|
case(LT_atom) |
136 |
|
|
m = 39.948E0_DP |
137 |
|
|
case(DZ_atom) |
138 |
|
|
m = 39.948E0_DP |
139 |
|
|
case default |
140 |
|
|
m = 0.0E0_DP |
141 |
|
|
end select |
142 |
|
|
|
143 |
|
|
end function atype_mass |
144 |
|
|
|
145 |
|
|
!integer returns -1 if error |
146 |
|
|
function atype_atoi(charident) result(t) |
147 |
|
|
|
148 |
|
|
character(len=2) :: charident |
149 |
|
|
integer :: t |
150 |
|
|
|
151 |
|
|
select case(charident) |
152 |
|
|
case('H') |
153 |
|
|
t = H_atom |
154 |
|
|
case('He') |
155 |
|
|
t = He_atom |
156 |
|
|
case('C') |
157 |
|
|
t = C_atom |
158 |
|
|
case('N') |
159 |
|
|
t = N_atom |
160 |
|
|
case('O') |
161 |
|
|
t = O_atom |
162 |
|
|
case('F') |
163 |
|
|
t = F_atom |
164 |
|
|
case('Ne') |
165 |
|
|
t = Ne_atom |
166 |
|
|
case('S') |
167 |
|
|
t = S_atom |
168 |
|
|
case('Cl') |
169 |
|
|
t = Cl_atom |
170 |
|
|
case('Ar') |
171 |
|
|
t = Ar_atom |
172 |
|
|
case('Ni') |
173 |
|
|
t = Ni_atom |
174 |
|
|
case('Cu') |
175 |
|
|
t = Cu_atom |
176 |
|
|
case('Rh') |
177 |
|
|
t = Rh_atom |
178 |
|
|
case('Pd') |
179 |
|
|
t = Pd_atom |
180 |
|
|
case('Ag') |
181 |
|
|
t = Ag_atom |
182 |
|
|
case('Ir') |
183 |
|
|
t = Ir_atom |
184 |
|
|
case('Pt') |
185 |
|
|
t = Pt_atom |
186 |
|
|
case('Au') |
187 |
|
|
t = Au_atom |
188 |
|
|
case('Pb') |
189 |
|
|
t = Pb_atom |
190 |
|
|
case('Al') |
191 |
|
|
t = Al_atom |
192 |
|
|
case('Br') |
193 |
|
|
t = Br_atom |
194 |
|
|
case('Kr') |
195 |
|
|
t = Kr_atom |
196 |
|
|
case('HV') |
197 |
|
|
t = HV_atom |
198 |
|
|
case('LT') |
199 |
|
|
t = LT_atom |
200 |
|
|
case('DZ') |
201 |
|
|
t = DZ_atom |
202 |
|
|
case default |
203 |
|
|
t = -1 |
204 |
|
|
end select |
205 |
|
|
|
206 |
|
|
end function atype_atoi |
207 |
|
|
|
208 |
|
|
! character returns XX if error. |
209 |
|
|
function atype_identifier(atype) result(c) |
210 |
|
|
|
211 |
|
|
integer :: atype |
212 |
|
|
character(len=2) :: c |
213 |
|
|
|
214 |
|
|
select case(atype) |
215 |
|
|
case(H_atom) |
216 |
|
|
c = 'H ' |
217 |
|
|
case(He_atom) |
218 |
|
|
c = 'He' |
219 |
|
|
case(C_atom) |
220 |
|
|
c = 'C ' |
221 |
|
|
case(N_atom) |
222 |
|
|
c = 'N ' |
223 |
|
|
case(O_atom) |
224 |
|
|
c = 'O ' |
225 |
|
|
case(F_atom) |
226 |
|
|
c = 'F ' |
227 |
|
|
case(Ne_atom) |
228 |
|
|
c = 'Ne' |
229 |
|
|
case(S_atom) |
230 |
|
|
c = 'S ' |
231 |
|
|
case(Cl_atom) |
232 |
|
|
c = 'Cl' |
233 |
|
|
case(Ar_atom) |
234 |
|
|
c = 'Ar' |
235 |
|
|
case(Ni_atom) |
236 |
|
|
c = 'Ni' |
237 |
|
|
case(Cu_atom) |
238 |
|
|
c = 'Cu' |
239 |
|
|
case(Rh_atom) |
240 |
|
|
c = 'Rh' |
241 |
|
|
case(Pd_atom) |
242 |
|
|
c = 'Pd' |
243 |
|
|
case(Ag_atom) |
244 |
|
|
c = 'Ag' |
245 |
|
|
case(Ir_atom) |
246 |
|
|
c = 'Ir' |
247 |
|
|
case(Pt_atom) |
248 |
|
|
c = 'Pt' |
249 |
|
|
case(Au_atom) |
250 |
|
|
c = 'Au' |
251 |
|
|
case(Pb_atom) |
252 |
|
|
c = 'Pb' |
253 |
|
|
case(Al_atom) |
254 |
|
|
c = 'Al' |
255 |
|
|
case(Br_atom) |
256 |
|
|
c = 'Br' |
257 |
|
|
case(Kr_atom) |
258 |
|
|
c = 'Kr' |
259 |
|
|
case(HV_atom) |
260 |
|
|
c = 'HV' |
261 |
|
|
case(LT_atom) |
262 |
|
|
c = 'LT' |
263 |
|
|
case(DZ_atom) |
264 |
|
|
c = 'DZ' |
265 |
|
|
case default |
266 |
|
|
c = 'XX' |
267 |
|
|
end select |
268 |
|
|
|
269 |
|
|
end function atype_identifier |
270 |
|
|
|
271 |
|
|
! integer |
272 |
|
|
function get_max_atype() result(max_atype) |
273 |
|
|
integer :: max_atype |
274 |
|
|
|
275 |
|
|
max_atype = natypes |
276 |
|
|
|
277 |
|
|
end function get_max_atype |
278 |
|
|
!-------------------------------------------------- |
279 |
|
|
!! multi atom models |
280 |
|
|
!-------------------------------------------------- |
281 |
|
|
function model_na(model) result(n) |
282 |
|
|
! returns the number of atoms in this type of molecule |
283 |
|
|
! assumes that any molecule it hasn't heard of before is simply an |
284 |
|
|
! atomic type (Ar, Au, etc.) (atomic types have only one atom per molecule) |
285 |
|
|
|
286 |
|
|
character(len=10) :: model |
287 |
|
|
integer :: n |
288 |
|
|
|
289 |
|
|
select case(model) |
290 |
|
|
case('o2') |
291 |
|
|
n = 2 |
292 |
|
|
case('n2') |
293 |
|
|
n = 2 |
294 |
|
|
case('br2') |
295 |
|
|
n = 2 |
296 |
|
|
case('cl2') |
297 |
|
|
n = 2 |
298 |
|
|
case('f2') |
299 |
|
|
n = 2 |
300 |
|
|
case('cs2') |
301 |
|
|
n = 3 |
302 |
|
|
case ('cfwater') |
303 |
|
|
n = 3 |
304 |
|
|
case default |
305 |
|
|
n = 1 |
306 |
|
|
end select |
307 |
|
|
end function model_na |
308 |
|
|
|
309 |
|
|
|
310 |
|
|
function model_atype(model, i) result(at) |
311 |
|
|
|
312 |
|
|
! returns the atom type of atom i in a molecule of type model |
313 |
|
|
! if it has never heard of this model type, and if i = 1, then |
314 |
|
|
! it assumes that model is the name of the atom, and uses the |
315 |
|
|
! atype functions to look up the atom type |
316 |
|
|
|
317 |
|
|
character(len=10) :: model |
318 |
|
|
character(len=10) :: tmp |
319 |
|
|
character(len=2) :: tmp2 |
320 |
|
|
integer i, at |
321 |
|
|
|
322 |
|
|
|
323 |
|
|
tmp = model |
324 |
|
|
tmp = tmp(1:machdep_lnblnk(tmp)) |
325 |
|
|
|
326 |
|
|
select case (tmp) |
327 |
|
|
case('h2') |
328 |
|
|
at = H_atom |
329 |
|
|
case('n2') |
330 |
|
|
at = N_atom |
331 |
|
|
case('o2') |
332 |
|
|
at = O_atom |
333 |
|
|
case('cl2') |
334 |
|
|
at = Cl_atom |
335 |
|
|
case('f2') |
336 |
|
|
at = F_atom |
337 |
|
|
case('br2') |
338 |
|
|
at = Br_atom |
339 |
|
|
case('cs2') |
340 |
|
|
if (i.eq.2) then |
341 |
|
|
at = C_atom |
342 |
|
|
else |
343 |
|
|
at = S_atom |
344 |
|
|
endif |
345 |
|
|
case('cfwater') |
346 |
|
|
if (i.eq.2) then |
347 |
|
|
at = O_atom |
348 |
|
|
else |
349 |
|
|
at = H_atom |
350 |
|
|
endif |
351 |
|
|
case default |
352 |
|
|
if (i.eq.1) then |
353 |
|
|
tmp2 = tmp(1:2) |
354 |
|
|
at = atype_atoi(tmp2) |
355 |
|
|
endif |
356 |
|
|
end select |
357 |
|
|
|
358 |
|
|
end function model_atype |
359 |
|
|
|
360 |
|
|
subroutine set_pos(molec, xcom, theta, phi, psi) |
361 |
|
|
use simulation |
362 |
|
|
integer :: j, i, molec, indx |
363 |
|
|
integer, dimension(3) :: lident |
364 |
|
|
|
365 |
|
|
real( kind = DP ), dimension(3) :: xcom, q_com |
366 |
|
|
real( kind = DP ) :: theta, phi, psi, pi |
367 |
|
|
real( kind = DP ), dimension(3) :: euler |
368 |
|
|
real( kind = DP ), dimension(3,3) :: rotate |
369 |
|
|
real( kind = DP ), dimension(3,3) :: q_standard |
370 |
|
|
real( kind = DP ), dimension(3) :: lmass |
371 |
|
|
real( kind = DP ) :: mtot |
372 |
|
|
real( kind = DP ), dimension(3,3) :: q_sites |
373 |
|
|
|
374 |
|
|
|
375 |
|
|
character(len=10) :: tmp |
376 |
|
|
|
377 |
|
|
|
378 |
|
|
pi = 4.0E0_DP * datan(1.0E0_DP) |
379 |
|
|
|
380 |
|
|
tmp = model(molec) |
381 |
|
|
tmp = tmp(1:machdep_lnblnk(tmp)) |
382 |
|
|
select case (tmp) |
383 |
|
|
case ('wca') |
384 |
|
|
q_standard(1,1) = 0.0E0_DP |
385 |
|
|
q_standard(2,1) = 0.0E0_DP |
386 |
|
|
q_standard(3,1) = 0.0E0_DP |
387 |
|
|
lmass(1) = 39.948E0_DP |
388 |
|
|
lident(1) = Ar_atom |
389 |
|
|
case ('cfwater') |
390 |
|
|
! Water: |
391 |
|
|
!---- O position |
392 |
|
|
q_standard(1,2) = -0.109554390183995380d-09 |
393 |
|
|
q_standard(2,2) = 0.208791443574830149E0_DP |
394 |
|
|
q_standard(3,2) = 0.000000000000000000d+00 |
395 |
|
|
lmass(2) = 15.9994E0_DP |
396 |
|
|
lident(2) = O_atom |
397 |
|
|
!---- H positions |
398 |
|
|
q_standard(1,1) = 0.757547679016284281E0_DP |
399 |
|
|
q_standard(2,1) = -0.592938967576377762E0_DP |
400 |
|
|
q_standard(3,1) = 0.00E0_DP |
401 |
|
|
lmass(1) = 1.0079 |
402 |
|
|
lident(1) = H_atom |
403 |
|
|
q_standard(1,3) = -0.757547678906731359E0_DP |
404 |
|
|
q_standard(2,3) = -0.592938967881253332E0_DP |
405 |
|
|
q_standard(3,3) = 0.00E0_DP |
406 |
|
|
lmass(3) = 1.0079 |
407 |
|
|
lident(3) = H_atom |
408 |
|
|
!---- End Water |
409 |
|
|
case ('cs2') |
410 |
|
|
! CS2: |
411 |
|
|
!---- C position |
412 |
|
|
q_standard(1,2) = 0.0E0_DP |
413 |
|
|
q_standard(2,2) = 0.0E0_DP |
414 |
|
|
q_standard(3,2) = 0.0E0_DP |
415 |
|
|
lmass(2) = 32.066E0_DP |
416 |
|
|
lident(2) = C_atom |
417 |
|
|
!---- S positions |
418 |
|
|
q_standard(1,1) = 1.57E0_DP |
419 |
|
|
q_standard(2,1) = 0.0E0_DP |
420 |
|
|
q_standard(3,1) = 0.00E0_DP |
421 |
|
|
lmass(1) = 12.011E0_DP |
422 |
|
|
lident(1) = S_atom |
423 |
|
|
q_standard(1,3) = -1.57E0_DP |
424 |
|
|
q_standard(2,3) = 0.0E0_DP |
425 |
|
|
q_standard(3,3) = 0.00E0_DP |
426 |
|
|
lmass(3) = 12.011E0_DP |
427 |
|
|
lident(3) = S_atom |
428 |
|
|
!---- End CS2 |
429 |
|
|
case ('cs2_fcc') |
430 |
|
|
! CS2: |
431 |
|
|
!---- C position |
432 |
|
|
q_standard(1,2) = 0.0E0_DP |
433 |
|
|
q_standard(2,2) = 0.0E0_DP |
434 |
|
|
q_standard(3,2) = 0.0E0_DP |
435 |
|
|
lmass(2) = 32.066E0_DP |
436 |
|
|
lident(2) = C_atom |
437 |
|
|
!---- S positions |
438 |
|
|
q_standard(1,1) = 1.57E0_DP |
439 |
|
|
q_standard(2,1) = 0.0E0_DP |
440 |
|
|
q_standard(3,1) = 0.00E0_DP |
441 |
|
|
lmass(1) = 12.011E0_DP |
442 |
|
|
lident(1) = S_atom |
443 |
|
|
q_standard(1,3) = -1.57E0_DP |
444 |
|
|
q_standard(2,3) = 0.0E0_DP |
445 |
|
|
q_standard(3,3) = 0.00E0_DP |
446 |
|
|
lmass(3) = 12.011E0_DP |
447 |
|
|
lident(3) = S_atom |
448 |
|
|
!---- End CS2 |
449 |
|
|
case default |
450 |
|
|
lident(1) = atype_atoi(tmp) |
451 |
|
|
lmass(1) = atype_mass(lident(1)) |
452 |
|
|
q_standard(1,1) = 0.0E0_DP |
453 |
|
|
q_standard(2,1) = 0.0E0_DP |
454 |
|
|
q_standard(3,1) = 0.0E0_DP |
455 |
|
|
end select |
456 |
|
|
|
457 |
|
|
! now adjust for COM: |
458 |
|
|
|
459 |
|
|
do j = 1, 3 |
460 |
|
|
q_com(j) = 0.0E0_DP |
461 |
|
|
enddo |
462 |
|
|
mtot = 0.0E0_DP |
463 |
|
|
do i = 1, na(molec) |
464 |
|
|
mtot = mtot + lmass(i) |
465 |
|
|
do j = 1, 3 |
466 |
|
|
q_com(j) = q_com(j) + lmass(i)*q_standard(j,i) |
467 |
|
|
enddo |
468 |
|
|
enddo |
469 |
|
|
|
470 |
|
|
do j = 1, 3 |
471 |
|
|
q_com(j) = q_com(j) / mtot |
472 |
|
|
enddo |
473 |
|
|
|
474 |
|
|
euler(1) = theta |
475 |
|
|
euler(2) = psi |
476 |
|
|
euler(3) = phi |
477 |
|
|
|
478 |
|
|
do i = 1, na(molec) |
479 |
|
|
do j = 1, 3 |
480 |
|
|
q_standard(j,i) = q_standard(j,i) - q_com(j) |
481 |
|
|
enddo |
482 |
|
|
enddo |
483 |
|
|
|
484 |
|
|
if (na(molec).ne.1) then |
485 |
|
|
call get_lab_pos(molec,xcom,euler,q_sites,q_standard,rotate) |
486 |
|
|
endif |
487 |
|
|
|
488 |
|
|
do i = 1, na(molec) |
489 |
|
|
indx = atom_index(molec,i) |
490 |
|
|
do j = 1, 3 |
491 |
|
|
mass(indx) = lmass(i) |
492 |
|
|
ident(indx) = lident(i) |
493 |
|
|
if (na(molec).eq.1) then |
494 |
|
|
q(j,indx) = q_standard(j,1) + xcom(j) |
495 |
|
|
else |
496 |
|
|
q(j,indx) = q_sites(j,i) |
497 |
|
|
endif |
498 |
|
|
enddo |
499 |
|
|
enddo |
500 |
|
|
|
501 |
|
|
return |
502 |
|
|
end subroutine set_pos |
503 |
|
|
|
504 |
|
|
|
505 |
|
|
subroutine get_lab_pos(molec, x, euler, q_sites, q_standard, & |
506 |
|
|
rotate) |
507 |
|
|
use simulation, ONLY: na |
508 |
|
|
real( kind = DP ), dimension(3) :: x,euler |
509 |
|
|
integer :: molec, j_site, i_dim |
510 |
|
|
|
511 |
|
|
! q and euler are the position and euler angles of the molecule |
512 |
|
|
|
513 |
|
|
real( kind = DP ), dimension(3,3) :: q_sites |
514 |
|
|
|
515 |
|
|
! q_sites(i,j) are the positions in the lab frame of the j-th |
516 |
|
|
! site on the molecule |
517 |
|
|
|
518 |
|
|
real( kind = DP ), dimension(3,3) :: q_standard |
519 |
|
|
|
520 |
|
|
! q_standard(i,j) is the unrotated position of the j-th site on a |
521 |
|
|
! standard water molecule |
522 |
|
|
|
523 |
|
|
real( kind = DP ), dimension(3,3) :: rotate |
524 |
|
|
real( kind = DP ), dimension(3) :: r_site_standard |
525 |
|
|
real( kind = DP ), dimension(3) :: r_site |
526 |
|
|
|
527 |
|
|
call deatom(euler(1),euler(2),euler(3),rotate) |
528 |
|
|
|
529 |
|
|
do j_site = 1, na(molec) |
530 |
|
|
do i_dim = 1, 3 |
531 |
|
|
r_site_standard(i_dim) = q_standard(i_dim,j_site) |
532 |
|
|
end do |
533 |
|
|
|
534 |
|
|
call vec_by_mat(r_site, rotate, r_site_standard) |
535 |
|
|
|
536 |
|
|
do i_dim = 1, 3 |
537 |
|
|
q_sites(i_dim, j_site) = x(i_dim) - r_site(i_dim) |
538 |
|
|
end do |
539 |
|
|
end do |
540 |
|
|
|
541 |
|
|
return |
542 |
|
|
end subroutine get_lab_pos |
543 |
|
|
|
544 |
|
|
|
545 |
|
|
subroutine vec_by_mat(vec_o,mat,vec_i) |
546 |
|
|
|
547 |
|
|
!---- vec_o := mat * vec_i |
548 |
|
|
|
549 |
|
|
integer :: i_dim,j_dim |
550 |
|
|
real( kind = DP ), dimension(3) :: vec_o |
551 |
|
|
real( kind = DP ), dimension(3,3) :: mat |
552 |
|
|
real( kind = DP ), dimension(3) :: vec_i |
553 |
|
|
|
554 |
|
|
do i_dim = 1,3 |
555 |
|
|
vec_o(i_dim) = 0.0E0_DP |
556 |
|
|
end do |
557 |
|
|
|
558 |
|
|
do i_dim = 1,3 |
559 |
|
|
do j_dim = 1,3 |
560 |
|
|
vec_o(j_dim) = vec_o(j_dim) + mat(j_dim,i_dim)*vec_i(i_dim) |
561 |
|
|
end do |
562 |
|
|
end do |
563 |
|
|
|
564 |
|
|
end subroutine vec_by_mat |
565 |
|
|
|
566 |
|
|
|
567 |
|
|
! H+ |
568 |
|
|
! Title : deatom_.f |
569 |
|
|
! Author : George Kaplan, from MIT matvec library |
570 |
|
|
! Date : 01 Mar 1990 |
571 |
|
|
! Synopsis: Generate rotation matrix from Euler angles |
572 |
|
|
! Keywords: ees, matvec, eatom, euler, matrix |
573 |
|
|
! @(#)deatom_.f 1.2 3/2/90 UCB SSL |
574 |
|
|
! Revisions: |
575 |
|
|
! mm/dd/yy name description |
576 |
|
|
! H- |
577 |
|
|
! U+ |
578 |
|
|
! Usage : call eatom(a, b, c, d) |
579 |
|
|
! Input : a, b, c: Euler angles |
580 |
|
|
! Output : d: Rotation matrix |
581 |
|
|
! U- |
582 |
|
|
! D+ |
583 |
|
|
! SUBROUTINE TO GENERATE ROTATION MATRIX GIVEN 3 EULER ANGLES |
584 |
|
|
! CALLING SEQUENCE |
585 |
|
|
! CALL EATOM(A,B,C,D) |
586 |
|
|
! (EATOM = EULER ANGLES TO MATRIX) |
587 |
|
|
! A FIRST EULER ANGLE |
588 |
|
|
! B SECOND EULER ANGLE |
589 |
|
|
! C THIRD EULER ANGLE |
590 |
|
|
! D RETURNED ROTATION MATRIX |
591 |
|
|
! |
592 |
|
|
! A,B,C ARE UNMODIFIED |
593 |
|
|
! |
594 |
|
|
! THE EULER ROTATION IS DEFINED AS FOLLOWS: |
595 |
|
|
! 1) ROTATE CCW BY ANGLE A ABOUT Z-AXIS GIVING X', Y', Z'-AXES |
596 |
|
|
! 2) ROTATE CCW BY ANGLE B ABOUT X'-AXIS GIVING X", Y", Z"-AXES |
597 |
|
|
! 3) ROTATE CCW BY ANGLE C ABOUT Z"-AXIS |
598 |
|
|
! |
599 |
|
|
! REFERENCE: |
600 |
|
|
! CLASSICAL MECHANICS, H. GOLDSTEIN, ADDISON WESLEY, 1950, P 107 FF |
601 |
|
|
! The three euler angles define a coordinate system rotation |
602 |
|
|
! from the original axes to a new set of axes. The rotation |
603 |
|
|
! matrix returned by EATOM rotates the coordinates of a vector |
604 |
|
|
! in the original axes to the coordinates in the new set. |
605 |
|
|
! x' = Dx |
606 |
|
|
! D- |
607 |
|
|
|
608 |
|
|
! ************************************************************ |
609 |
|
|
|
610 |
|
|
SUBROUTINE deatom(A,B,C,D) |
611 |
|
|
real( kind = DP ) :: a,b,c |
612 |
|
|
real( kind = DP ), dimension(9) :: d |
613 |
|
|
real( kind = DP ) :: sina, cosa, sinb, cosb, sinc, cosc |
614 |
|
|
real( kind = DP ) :: exp1, exp2 |
615 |
|
|
|
616 |
|
|
SINA=sin(A) |
617 |
|
|
COSA=cos(A) |
618 |
|
|
|
619 |
|
|
SINB=sin(B) |
620 |
|
|
COSB=cos(B) |
621 |
|
|
|
622 |
|
|
SINC=sin(C) |
623 |
|
|
COSC=cos(C) |
624 |
|
|
|
625 |
|
|
EXP1=SINA*COSB |
626 |
|
|
EXP2=COSA*COSB |
627 |
|
|
|
628 |
|
|
D(1)=+COSA*COSC-EXP1*SINC |
629 |
|
|
D(2)=-COSA*SINC-EXP1*COSC |
630 |
|
|
D(3)=+SINA*SINB |
631 |
|
|
|
632 |
|
|
D(4)=+SINA*COSC+EXP2*SINC |
633 |
|
|
D(5)=-SINA*SINC+EXP2*COSC |
634 |
|
|
D(6)=-COSA*SINB |
635 |
|
|
|
636 |
|
|
D(7)=SINB*SINC |
637 |
|
|
D(8)=SINB*COSC |
638 |
|
|
D(9)=COSB |
639 |
|
|
|
640 |
|
|
END SUBROUTINE deatom |
641 |
|
|
|
642 |
|
|
|
643 |
|
|
|
644 |
|
|
end module model_module |