1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> Torsions added (DG) 05-Dec-2000 9!> Variable names changed (DG) 05-Dec-2000 10!> \author CJM 11! ************************************************************************************************** 12MODULE fist_intra_force 13 14 USE atprop_types, ONLY: atprop_type 15 USE cell_types, ONLY: cell_type,& 16 pbc 17 USE cp_log_handling, ONLY: cp_get_default_logger,& 18 cp_logger_type 19 USE distribution_1d_types, ONLY: distribution_1d_type 20 USE kinds, ONLY: dp 21 USE mol_force, ONLY: force_bends,& 22 force_bonds,& 23 force_imp_torsions,& 24 force_opbends,& 25 force_torsions,& 26 get_pv_bend,& 27 get_pv_bond,& 28 get_pv_torsion 29 USE molecule_kind_types, ONLY: & 30 bend_type, bond_type, get_molecule_kind, impr_type, molecule_kind_type, opbend_type, & 31 shell_type, torsion_type, ub_type 32 USE molecule_types, ONLY: get_molecule,& 33 molecule_type 34 USE particle_types, ONLY: particle_type 35#include "./base/base_uses.f90" 36 37 IMPLICIT NONE 38 39 PRIVATE 40 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_intra_force' 41 PUBLIC :: force_intra_control 42 43CONTAINS 44 45! ************************************************************************************************** 46!> \brief Computes the the intramolecular energies, forces, and pressure tensors 47!> \param molecule_set ... 48!> \param molecule_kind_set ... 49!> \param local_molecules ... 50!> \param particle_set ... 51!> \param shell_particle_set ... 52!> \param core_particle_set ... 53!> \param pot_bond ... 54!> \param pot_bend ... 55!> \param pot_urey_bradley ... 56!> \param pot_torsion ... 57!> \param pot_imp_torsion ... 58!> \param pot_opbend ... 59!> \param pot_shell ... 60!> \param pv_bond ... 61!> \param pv_bend ... 62!> \param pv_urey_bradley ... 63!> \param pv_torsion ... 64!> \param pv_imp_torsion ... 65!> \param pv_opbend ... 66!> \param f_bond ... 67!> \param f_bend ... 68!> \param f_torsion ... 69!> \param f_ub ... 70!> \param f_imptor ... 71!> \param f_opbend ... 72!> \param cell ... 73!> \param use_virial ... 74!> \param atprop_env ... 75!> \par History 76!> none 77!> \author CJM 78! ************************************************************************************************** 79 SUBROUTINE force_intra_control(molecule_set, molecule_kind_set, & 80 local_molecules, particle_set, shell_particle_set, core_particle_set, & 81 pot_bond, pot_bend, pot_urey_bradley, pot_torsion, pot_imp_torsion, & 82 pot_opbend, pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, & 83 pv_imp_torsion, pv_opbend, f_bond, f_bend, f_torsion, f_ub, & 84 f_imptor, f_opbend, cell, use_virial, atprop_env) 85 86 TYPE(molecule_type), POINTER :: molecule_set(:) 87 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 88 TYPE(distribution_1d_type), POINTER :: local_molecules 89 TYPE(particle_type), POINTER :: particle_set(:), shell_particle_set(:), & 90 core_particle_set(:) 91 REAL(KIND=dp), INTENT(INOUT) :: pot_bond, pot_bend, pot_urey_bradley, & 92 pot_torsion, pot_imp_torsion, & 93 pot_opbend, pot_shell 94 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pv_bond, pv_bend, pv_urey_bradley, & 95 pv_torsion, pv_imp_torsion, pv_opbend 96 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 97 OPTIONAL :: f_bond, f_bend, f_torsion, f_ub, & 98 f_imptor, f_opbend 99 TYPE(cell_type), POINTER :: cell 100 LOGICAL, INTENT(IN) :: use_virial 101 TYPE(atprop_type), POINTER :: atprop_env 102 103 CHARACTER(len=*), PARAMETER :: routineN = 'force_intra_control', & 104 routineP = moduleN//':'//routineN 105 106 INTEGER :: first_atom, handle, i, ibend, ibond, ikind, imol, imul, index_a, index_b, & 107 index_c, index_d, iopbend, ishell, itorsion, nbends, nbonds, nimptors, nkind, & 108 nmol_per_kind, nopbends, nshell, ntorsions, nub 109 LOGICAL :: atener, atstress 110 REAL(KIND=dp) :: d12, d32, dist, dist1, dist2, energy, & 111 fscalar, id12, id32, is32, ism, isn, & 112 k2_spring, k4_spring, r2, s32, sm, sn, & 113 theta 114 REAL(KIND=dp), DIMENSION(3) :: b12, b32, g1, g2, g3, gt1, gt2, gt3, & 115 gt4, k1, k2, k3, k4, rij, t12, t32, & 116 t34, t41, t42, t43, tm, tn 117 REAL(KIND=dp), DIMENSION(:), POINTER :: ener_a 118 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pv_a 119 TYPE(bend_type), POINTER :: bend_list(:) 120 TYPE(bond_type), POINTER :: bond_list(:) 121 TYPE(cp_logger_type), POINTER :: logger 122 TYPE(impr_type), POINTER :: impr_list(:) 123 TYPE(molecule_kind_type), POINTER :: molecule_kind 124 TYPE(molecule_type), POINTER :: molecule 125 TYPE(opbend_type), POINTER :: opbend_list(:) 126 TYPE(shell_type), POINTER :: shell_list(:) 127 TYPE(torsion_type), POINTER :: torsion_list(:) 128 TYPE(ub_type), POINTER :: ub_list(:) 129 130 CALL timeset(routineN, handle) 131 NULLIFY (logger) 132 logger => cp_get_default_logger() 133 134 IF (PRESENT(f_bond)) f_bond = 0.0_dp 135 IF (PRESENT(f_bend)) f_bend = 0.0_dp 136 IF (PRESENT(f_torsion)) f_torsion = 0.0_dp 137 IF (PRESENT(f_imptor)) f_imptor = 0.0_dp 138 IF (PRESENT(f_ub)) f_ub = 0.0_dp 139 140 pot_bond = 0.0_dp 141 pot_bend = 0.0_dp 142 pot_urey_bradley = 0.0_dp 143 pot_torsion = 0.0_dp 144 pot_imp_torsion = 0.0_dp 145 pot_opbend = 0.0_dp 146 pot_shell = 0.0_dp 147 148 atener = atprop_env%energy 149 IF (atener) ener_a => atprop_env%atener 150 atstress = atprop_env%stress 151 IF (atstress) pv_a => atprop_env%atstress 152 153 nkind = SIZE(molecule_kind_set) 154 MOL: DO ikind = 1, nkind 155 nmol_per_kind = local_molecules%n_el(ikind) 156 157 DO imol = 1, nmol_per_kind 158 i = local_molecules%list(ikind)%array(imol) 159 molecule => molecule_set(i) 160 molecule_kind => molecule%molecule_kind 161 CALL get_molecule_kind(molecule_kind, nbend=nbends, nbond=nbonds, & 162 nimpr=nimptors, nub=nub, ntorsion=ntorsions, & 163 nopbend=nopbends, nshell=nshell, & 164 bond_list=bond_list, ub_list=ub_list, & 165 bend_list=bend_list, torsion_list=torsion_list, & 166 impr_list=impr_list, opbend_list=opbend_list, & 167 shell_list=shell_list) 168 169 CALL get_molecule(molecule, first_atom=first_atom) 170 171 BOND: DO ibond = 1, nbonds 172 index_a = bond_list(ibond)%a + first_atom - 1 173 index_b = bond_list(ibond)%b + first_atom - 1 174 rij = particle_set(index_a)%r - particle_set(index_b)%r 175 rij = pbc(rij, cell) 176 CALL force_bonds(bond_list(ibond)%bond_kind%id_type, rij, & 177 bond_list(ibond)%bond_kind%r0, & 178 bond_list(ibond)%bond_kind%k, & 179 bond_list(ibond)%bond_kind%cs, & 180 energy, fscalar) 181 pot_bond = pot_bond + energy 182 IF (atener) THEN 183 ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy 184 ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy 185 END IF 186 187 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar 188 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar 189 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar 190 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar 191 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar 192 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar 193 194 ! computing the pressure tensor 195 k2 = -rij*fscalar 196 IF (use_virial) CALL get_pv_bond(k2, rij, pv_bond) 197 IF (atstress) THEN 198 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_a)) 199 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_b)) 200 END IF 201 202 ! the contribution from the bonds. ONLY FOR DEBUG 203 IF (PRESENT(f_bond)) THEN 204 f_bond(1, index_a) = f_bond(1, index_a) - rij(1)*fscalar 205 f_bond(2, index_a) = f_bond(2, index_a) - rij(2)*fscalar 206 f_bond(3, index_a) = f_bond(3, index_a) - rij(3)*fscalar 207 f_bond(1, index_b) = f_bond(1, index_b) + rij(1)*fscalar 208 f_bond(2, index_b) = f_bond(2, index_b) + rij(2)*fscalar 209 f_bond(3, index_b) = f_bond(3, index_b) + rij(3)*fscalar 210 END IF 211 212 END DO BOND 213 214 SHELL: DO ishell = 1, nshell 215 index_a = shell_list(ishell)%a + first_atom - 1 216 index_b = particle_set(index_a)%shell_index 217 rij = core_particle_set(index_b)%r - shell_particle_set(index_b)%r 218 rij = pbc(rij, cell) 219 k2_spring = shell_list(ishell)%shell_kind%k2_spring 220 k4_spring = shell_list(ishell)%shell_kind%k4_spring 221 r2 = DOT_PRODUCT(rij, rij) 222 energy = 0.5_dp*(k2_spring + k4_spring*r2/12.0_dp)*r2 223 fscalar = k2_spring + k4_spring*r2/6.0_dp 224 pot_shell = pot_shell + energy 225 IF (atener) THEN 226 ener_a(index_a) = ener_a(index_a) + energy 227 END IF 228 core_particle_set(index_b)%f(1) = core_particle_set(index_b)%f(1) - rij(1)*fscalar 229 core_particle_set(index_b)%f(2) = core_particle_set(index_b)%f(2) - rij(2)*fscalar 230 core_particle_set(index_b)%f(3) = core_particle_set(index_b)%f(3) - rij(3)*fscalar 231 shell_particle_set(index_b)%f(1) = shell_particle_set(index_b)%f(1) + rij(1)*fscalar 232 shell_particle_set(index_b)%f(2) = shell_particle_set(index_b)%f(2) + rij(2)*fscalar 233 shell_particle_set(index_b)%f(3) = shell_particle_set(index_b)%f(3) + rij(3)*fscalar 234 ! Compute the pressure tensor, if requested 235 IF (use_virial) THEN 236 k1 = -rij*fscalar 237 CALL get_pv_bond(k1, rij, pv_bond) 238 END IF 239 IF (atstress) THEN 240 CALL get_pv_bond(k1, rij, pv_a(:, :, index_a)) 241 END IF 242 END DO SHELL 243 244 UREY_BRADLEY: DO ibend = 1, nub 245 index_a = ub_list(ibend)%a + first_atom - 1 246 index_b = ub_list(ibend)%c + first_atom - 1 247 rij = particle_set(index_a)%r - particle_set(index_b)%r 248 rij = pbc(rij, cell) 249 CALL force_bonds(ub_list(ibend)%ub_kind%id_type, rij, & 250 ub_list(ibend)%ub_kind%r0, & 251 ub_list(ibend)%ub_kind%k, 0.0_dp, energy, fscalar) 252 pot_urey_bradley = pot_urey_bradley + energy 253 IF (atener) THEN 254 ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy 255 ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy 256 END IF 257 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar 258 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar 259 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar 260 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar 261 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar 262 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar 263 264 ! computing the pressure tensor 265 k2 = -rij*fscalar 266 IF (use_virial) CALL get_pv_bond(k2, rij, pv_urey_bradley) 267 IF (atstress) THEN 268 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_a)) 269 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_b)) 270 END IF 271 272 ! the contribution from the ub. ONLY FOR DEBUG 273 IF (PRESENT(f_ub)) THEN 274 f_ub(:, index_a) = f_ub(:, index_a) - rij*fscalar 275 f_ub(:, index_b) = f_ub(:, index_b) + rij*fscalar 276 END IF 277 278 END DO UREY_BRADLEY 279 280 BEND: DO ibend = 1, nbends 281 index_a = bend_list(ibend)%a + first_atom - 1 282 index_b = bend_list(ibend)%b + first_atom - 1 283 index_c = bend_list(ibend)%c + first_atom - 1 284 b12 = particle_set(index_a)%r - particle_set(index_b)%r 285 b32 = particle_set(index_c)%r - particle_set(index_b)%r 286 b12 = pbc(b12, cell) 287 b32 = pbc(b32, cell) 288 d12 = SQRT(DOT_PRODUCT(b12, b12)) 289 id12 = 1.0_dp/d12 290 d32 = SQRT(DOT_PRODUCT(b32, b32)) 291 id32 = 1.0_dp/d32 292 dist = DOT_PRODUCT(b12, b32) 293 theta = (dist*id12*id32) 294 IF (theta < -1.0_dp) theta = -1.0_dp 295 IF (theta > +1.0_dp) theta = +1.0_dp 296 theta = ACOS(theta) 297 CALL force_bends(bend_list(ibend)%bend_kind%id_type, & 298 b12, b32, d12, d32, id12, id32, dist, theta, & 299 bend_list(ibend)%bend_kind%theta0, & 300 bend_list(ibend)%bend_kind%k, & 301 bend_list(ibend)%bend_kind%cb, & 302 bend_list(ibend)%bend_kind%r012, & 303 bend_list(ibend)%bend_kind%r032, & 304 bend_list(ibend)%bend_kind%kbs12, & 305 bend_list(ibend)%bend_kind%kbs32, & 306 bend_list(ibend)%bend_kind%kss, & 307 bend_list(ibend)%bend_kind%legendre, & 308 g1, g2, g3, energy, fscalar) 309 pot_bend = pot_bend + energy 310 IF (atener) THEN 311 ener_a(index_a) = ener_a(index_a) + energy/3._dp 312 ener_a(index_b) = ener_a(index_b) + energy/3._dp 313 ener_a(index_c) = ener_a(index_c) + energy/3._dp 314 END IF 315 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + g1(1)*fscalar 316 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + g1(2)*fscalar 317 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + g1(3)*fscalar 318 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + g2(1)*fscalar 319 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + g2(2)*fscalar 320 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + g2(3)*fscalar 321 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + g3(1)*fscalar 322 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + g3(2)*fscalar 323 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + g3(3)*fscalar 324 325 ! computing the pressure tensor 326 k1 = fscalar*g1 327 k3 = fscalar*g3 328 IF (use_virial) CALL get_pv_bend(k1, k3, b12, b32, pv_bend) 329 IF (atstress) THEN 330 k1 = fscalar*g1/3._dp 331 k3 = fscalar*g3/3._dp 332 CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_a)) 333 CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_b)) 334 CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_c)) 335 END IF 336 337 ! the contribution from the bends. ONLY FOR DEBUG 338 IF (PRESENT(f_bend)) THEN 339 f_bend(:, index_a) = f_bend(:, index_a) + fscalar*g1 340 f_bend(:, index_b) = f_bend(:, index_b) + fscalar*g2 341 f_bend(:, index_c) = f_bend(:, index_c) + fscalar*g3 342 END IF 343 END DO BEND 344 345 TORSION: DO itorsion = 1, ntorsions 346 index_a = torsion_list(itorsion)%a + first_atom - 1 347 index_b = torsion_list(itorsion)%b + first_atom - 1 348 index_c = torsion_list(itorsion)%c + first_atom - 1 349 index_d = torsion_list(itorsion)%d + first_atom - 1 350 t12 = particle_set(index_a)%r - particle_set(index_b)%r 351 t32 = particle_set(index_c)%r - particle_set(index_b)%r 352 t34 = particle_set(index_c)%r - particle_set(index_d)%r 353 t43 = particle_set(index_d)%r - particle_set(index_c)%r 354 t12 = pbc(t12, cell) 355 t32 = pbc(t32, cell) 356 t34 = pbc(t34, cell) 357 t43 = pbc(t43, cell) 358 ! t12 x t32 359 tm(1) = t12(2)*t32(3) - t32(2)*t12(3) 360 tm(2) = -t12(1)*t32(3) + t32(1)*t12(3) 361 tm(3) = t12(1)*t32(2) - t32(1)*t12(2) 362 ! t32 x t34 363 tn(1) = t32(2)*t34(3) - t34(2)*t32(3) 364 tn(2) = -t32(1)*t34(3) + t34(1)*t32(3) 365 tn(3) = t32(1)*t34(2) - t34(1)*t32(2) 366 sm = SQRT(DOT_PRODUCT(tm, tm)) 367 ism = 1.0_dp/sm 368 sn = SQRT(DOT_PRODUCT(tn, tn)) 369 isn = 1.0_dp/sn 370 s32 = SQRT(DOT_PRODUCT(t32, t32)) 371 is32 = 1.0_dp/s32 372 dist1 = DOT_PRODUCT(t12, t32) 373 dist2 = DOT_PRODUCT(t34, t32) 374 DO imul = 1, torsion_list(itorsion)%torsion_kind%nmul 375 CALL force_torsions(torsion_list(itorsion)%torsion_kind%id_type, & 376 s32, is32, ism, isn, dist1, dist2, tm, tn, t12, & 377 torsion_list(itorsion)%torsion_kind%k(imul), & 378 torsion_list(itorsion)%torsion_kind%phi0(imul), & 379 torsion_list(itorsion)%torsion_kind%m(imul), & 380 gt1, gt2, gt3, gt4, energy, fscalar) 381 pot_torsion = pot_torsion + energy 382 IF (atener) THEN 383 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp 384 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp 385 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp 386 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp 387 END IF 388 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar 389 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar 390 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar 391 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar 392 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar 393 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar 394 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar 395 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar 396 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar 397 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar 398 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar 399 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar 400 401 ! computing the pressure tensor 402 k1 = fscalar*gt1 403 k3 = fscalar*gt3 404 k4 = fscalar*gt4 405 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_torsion) 406 IF (atstress) THEN 407 k1 = 0.25_dp*fscalar*gt1 408 k3 = 0.25_dp*fscalar*gt3 409 k4 = 0.25_dp*fscalar*gt4 410 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a)) 411 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b)) 412 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c)) 413 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_d)) 414 END IF 415 416 ! the contribution from the torsions. ONLY FOR DEBUG 417 IF (PRESENT(f_torsion)) THEN 418 f_torsion(:, index_a) = f_torsion(:, index_a) + fscalar*gt1 419 f_torsion(:, index_b) = f_torsion(:, index_b) + fscalar*gt2 420 f_torsion(:, index_c) = f_torsion(:, index_c) + fscalar*gt3 421 f_torsion(:, index_d) = f_torsion(:, index_d) + fscalar*gt4 422 END IF 423 END DO 424 END DO TORSION 425 426 IMP_TORSION: DO itorsion = 1, nimptors 427 index_a = impr_list(itorsion)%a + first_atom - 1 428 index_b = impr_list(itorsion)%b + first_atom - 1 429 index_c = impr_list(itorsion)%c + first_atom - 1 430 index_d = impr_list(itorsion)%d + first_atom - 1 431 t12 = particle_set(index_a)%r - particle_set(index_b)%r 432 t32 = particle_set(index_c)%r - particle_set(index_b)%r 433 t34 = particle_set(index_c)%r - particle_set(index_d)%r 434 t43 = particle_set(index_d)%r - particle_set(index_c)%r 435 t12 = pbc(t12, cell) 436 t32 = pbc(t32, cell) 437 t34 = pbc(t34, cell) 438 t43 = pbc(t43, cell) 439 ! t12 x t32 440 tm(1) = t12(2)*t32(3) - t32(2)*t12(3) 441 tm(2) = -t12(1)*t32(3) + t32(1)*t12(3) 442 tm(3) = t12(1)*t32(2) - t32(1)*t12(2) 443 ! t32 x t34 444 tn(1) = t32(2)*t34(3) - t34(2)*t32(3) 445 tn(2) = -t32(1)*t34(3) + t34(1)*t32(3) 446 tn(3) = t32(1)*t34(2) - t34(1)*t32(2) 447 sm = SQRT(DOT_PRODUCT(tm, tm)) 448 ism = 1.0_dp/sm 449 sn = SQRT(DOT_PRODUCT(tn, tn)) 450 isn = 1.0_dp/sn 451 s32 = SQRT(DOT_PRODUCT(t32, t32)) 452 is32 = 1.0_dp/s32 453 dist1 = DOT_PRODUCT(t12, t32) 454 dist2 = DOT_PRODUCT(t34, t32) 455 CALL force_imp_torsions(impr_list(itorsion)%impr_kind%id_type, & 456 s32, is32, ism, isn, dist1, dist2, tm, tn, t12, & 457 impr_list(itorsion)%impr_kind%k, & 458 impr_list(itorsion)%impr_kind%phi0, & 459 gt1, gt2, gt3, gt4, energy, fscalar) 460 pot_imp_torsion = pot_imp_torsion + energy 461 IF (atener) THEN 462 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp 463 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp 464 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp 465 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp 466 END IF 467 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar 468 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar 469 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar 470 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar 471 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar 472 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar 473 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar 474 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar 475 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar 476 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar 477 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar 478 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar 479 480 ! computing the pressure tensor 481 k1 = fscalar*gt1 482 k3 = fscalar*gt3 483 k4 = fscalar*gt4 484 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_imp_torsion) 485 IF (atstress) THEN 486 k1 = 0.25_dp*fscalar*gt1 487 k3 = 0.25_dp*fscalar*gt3 488 k4 = 0.25_dp*fscalar*gt4 489 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a)) 490 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b)) 491 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c)) 492 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_d)) 493 END IF 494 495 ! the contribution from the torsions. ONLY FOR DEBUG 496 IF (PRESENT(f_imptor)) THEN 497 f_imptor(:, index_a) = f_imptor(:, index_a) + fscalar*gt1 498 f_imptor(:, index_b) = f_imptor(:, index_b) + fscalar*gt2 499 f_imptor(:, index_c) = f_imptor(:, index_c) + fscalar*gt3 500 f_imptor(:, index_d) = f_imptor(:, index_d) + fscalar*gt4 501 END IF 502 END DO IMP_TORSION 503 504 OPBEND: DO iopbend = 1, nopbends 505 index_a = opbend_list(iopbend)%a + first_atom - 1 506 index_b = opbend_list(iopbend)%b + first_atom - 1 507 index_c = opbend_list(iopbend)%c + first_atom - 1 508 index_d = opbend_list(iopbend)%d + first_atom - 1 509 510 t12 = particle_set(index_a)%r - particle_set(index_b)%r 511 t32 = particle_set(index_c)%r - particle_set(index_b)%r 512 t34 = particle_set(index_c)%r - particle_set(index_d)%r 513 t43 = particle_set(index_d)%r - particle_set(index_c)%r 514 t41 = particle_set(index_d)%r - particle_set(index_a)%r 515 t42 = pbc(t41 + t12, cell) 516 t12 = pbc(t12, cell) 517 t32 = pbc(t32, cell) 518 t41 = pbc(t41, cell) 519 t43 = pbc(t43, cell) 520 ! tm = t32 x t12 521 tm(1) = t32(2)*t12(3) - t12(2)*t32(3) 522 tm(2) = -t32(1)*t12(3) + t12(1)*t32(3) 523 tm(3) = t32(1)*t12(2) - t12(1)*t32(2) 524 sm = SQRT(DOT_PRODUCT(tm, tm)) 525 s32 = SQRT(DOT_PRODUCT(t32, t32)) 526 CALL force_opbends(opbend_list(iopbend)%opbend_kind%id_type, & 527 s32, tm, t41, t42, t43, & 528 opbend_list(iopbend)%opbend_kind%k, & 529 opbend_list(iopbend)%opbend_kind%phi0, & 530 gt1, gt2, gt3, gt4, energy, fscalar) 531 pot_opbend = pot_opbend + energy 532 IF (atener) THEN 533 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp 534 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp 535 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp 536 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp 537 END IF 538 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar 539 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar 540 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar 541 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar 542 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar 543 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar 544 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar 545 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar 546 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar 547 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar 548 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar 549 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar 550 551 ! computing the pressure tensor 552 k1 = fscalar*gt1 553 k3 = fscalar*gt3 554 k4 = fscalar*gt4 555 556 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_opbend) 557 IF (atstress) THEN 558 k1 = 0.25_dp*fscalar*gt1 559 k3 = 0.25_dp*fscalar*gt3 560 k4 = 0.25_dp*fscalar*gt4 561 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a)) 562 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b)) 563 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c)) 564 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_d)) 565 END IF 566 567 ! the contribution from the opbends. ONLY FOR DEBUG 568 IF (PRESENT(f_opbend)) THEN 569 f_opbend(:, index_a) = f_opbend(:, index_a) + fscalar*gt1 570 f_opbend(:, index_b) = f_opbend(:, index_b) + fscalar*gt2 571 f_opbend(:, index_c) = f_opbend(:, index_c) + fscalar*gt3 572 f_opbend(:, index_d) = f_opbend(:, index_d) + fscalar*gt4 573 END IF 574 END DO OPBEND 575 END DO 576 END DO MOL 577 578 CALL timestop(handle) 579 580 END SUBROUTINE force_intra_control 581 582END MODULE fist_intra_force 583 584