1 |
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 |