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