1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> Torsions added (DG) 05-Dec-2000 9!> \author CJM 10! ************************************************************************************************** 11MODULE mol_force 12 13 USE force_field_kind_types, ONLY: & 14 do_ff_amber, do_ff_charmm, do_ff_cubic, do_ff_fues, do_ff_g87, do_ff_g96, do_ff_harmonic, & 15 do_ff_legendre, do_ff_mixed_bend_stretch, do_ff_mm2, do_ff_mm3, do_ff_mm4, do_ff_morse, & 16 do_ff_opls, do_ff_quartic, legendre_data_type 17 USE kinds, ONLY: dp 18#include "./base/base_uses.f90" 19 20 IMPLICIT NONE 21 22 PRIVATE 23 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mol_force' 24 PUBLIC :: force_bonds, force_bends, force_torsions, force_imp_torsions, force_opbends 25 PUBLIC :: get_pv_bond, get_pv_bend, get_pv_torsion 26 27CONTAINS 28 29! ************************************************************************************************** 30!> \brief Computes the forces from the bonds 31!> \param id_type ... 32!> \param rij ... 33!> \param r0 ... 34!> \param k ... 35!> \param cs ... 36!> \param energy ... 37!> \param fscalar ... 38!> \author CJM 39! ************************************************************************************************** 40 SUBROUTINE force_bonds(id_type, rij, r0, k, cs, energy, fscalar) 41 INTEGER, INTENT(IN) :: id_type 42 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rij 43 REAL(KIND=dp), INTENT(IN) :: r0, k(3), cs 44 REAL(KIND=dp), INTENT(OUT) :: energy, fscalar 45 46 CHARACTER(len=*), PARAMETER :: routineN = 'force_bonds', routineP = moduleN//':'//routineN 47 REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp, & 48 f13 = 1.0_dp/3.0_dp, & 49 f14 = 1.0_dp/4.0_dp 50 51 REAL(KIND=dp) :: dij, disp 52 53 SELECT CASE (id_type) 54 CASE (do_ff_quartic) 55 dij = SQRT(DOT_PRODUCT(rij, rij)) 56 disp = dij - r0 57 energy = (f12*k(1) + (f13*k(2) + f14*k(3)*disp)*disp)*disp*disp 58 fscalar = ((k(1) + (k(2) + k(3)*disp)*disp)*disp)/dij 59 CASE (do_ff_morse) 60 dij = SQRT(DOT_PRODUCT(rij, rij)) 61 disp = dij - r0 62 energy = k(1)*((1 - EXP(-k(2)*disp))**2 - 1) 63 fscalar = 2*k(1)*k(2)*EXP(-k(2)*disp)*(1 - EXP(-k(2)*disp))/dij 64 CASE (do_ff_cubic) 65 dij = SQRT(DOT_PRODUCT(rij, rij)) 66 disp = dij - r0 67 energy = k(1)*disp**2*(1 + cs*disp + 7.0_dp/12.0_dp*cs**2*disp**2) 68 fscalar = (2.0_dp*k(1)*disp*(1 + cs*disp + 7.0_dp/12.0_dp*cs**2*disp**2) + & 69 k(1)*disp**2*(cs + 2.0_dp*7.0_dp/12.0_dp*cs**2*disp))/dij 70 CASE (do_ff_g96) 71 ! From GROMOS... 72 ! V = (1/4)*Kb*(rij**2 - bij**2)**2 73 dij = DOT_PRODUCT(rij, rij) 74 disp = dij - r0*r0 75 energy = f14*k(1)*disp*disp 76 fscalar = k(1)*disp 77 CASE (do_ff_charmm, do_ff_amber) 78 dij = SQRT(DOT_PRODUCT(rij, rij)) 79 disp = dij - r0 80 IF (ABS(disp) < EPSILON(1.0_dp)) THEN 81 energy = 0.0_dp 82 fscalar = 0.0_dp 83 ELSE 84 energy = k(1)*disp*disp 85 fscalar = 2.0_dp*k(1)*disp/dij 86 END IF 87 CASE (do_ff_harmonic, do_ff_g87) 88 dij = SQRT(DOT_PRODUCT(rij, rij)) 89 disp = dij - r0 90 IF (ABS(disp) < EPSILON(1.0_dp)) THEN 91 energy = 0.0_dp 92 fscalar = 0.0_dp 93 ELSE 94 energy = f12*k(1)*disp*disp 95 fscalar = k(1)*disp/dij 96 END IF 97 CASE (do_ff_fues) 98 dij = SQRT(DOT_PRODUCT(rij, rij)) 99 disp = r0/dij 100 energy = f12*k(1)*r0*r0*(1.0_dp + disp*(disp - 2.0_dp)) 101 fscalar = k(1)*r0*disp*disp*(1.0_dp - disp)/dij 102 CASE DEFAULT 103 CPABORT("Unmatched bond kind") 104 END SELECT 105 106 END SUBROUTINE force_bonds 107 108! ************************************************************************************************** 109!> \brief Computes the forces from the bends 110!> \param id_type ... 111!> \param b12 ... 112!> \param b32 ... 113!> \param d12 ... 114!> \param d32 ... 115!> \param id12 ... 116!> \param id32 ... 117!> \param dist ... 118!> \param theta ... 119!> \param theta0 ... 120!> \param k ... 121!> \param cb ... 122!> \param r012 ... 123!> \param r032 ... 124!> \param kbs12 ... 125!> \param kbs32 ... 126!> \param kss ... 127!> \param legendre ... 128!> \param g1 ... 129!> \param g2 ... 130!> \param g3 ... 131!> \param energy ... 132!> \param fscalar ... 133!> \par History 134!> Legendre Bonding Potential added 2015-11-02 [Felix Uhl] 135!> \author CJM 136! ************************************************************************************************** 137 SUBROUTINE force_bends(id_type, b12, b32, d12, d32, id12, id32, dist, & 138 theta, theta0, k, cb, r012, r032, kbs12, kbs32, kss, legendre, g1, g2, g3, energy, fscalar) 139 INTEGER, INTENT(IN) :: id_type 140 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: b12, b32 141 REAL(KIND=dp), INTENT(IN) :: d12, d32, id12, id32, dist, theta, & 142 theta0, k, cb, r012, r032, kbs12, & 143 kbs32, kss 144 TYPE(legendre_data_type), INTENT(IN) :: legendre 145 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: g1, g2, g3 146 REAL(KIND=dp), INTENT(OUT) :: energy, fscalar 147 148 CHARACTER(len=*), PARAMETER :: routineN = 'force_bends', routineP = moduleN//':'//routineN 149 REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp 150 151 INTEGER :: i 152 REAL(KIND=dp) :: ctheta, denom, disp12, disp32, y0, y1, & 153 y2, yd0, yd1, yd2 154 155 SELECT CASE (id_type) 156 CASE (do_ff_g96) 157 energy = f12*k*(COS(theta) - theta0)**2 158 fscalar = -k*(COS(theta) - theta0) 159 g1 = (b32*id32 - b12*id12*COS(theta))*id12 160 g3 = (b12*id12 - b32*id32*COS(theta))*id32 161 g2 = -g1 - g3 162 CASE (do_ff_charmm, do_ff_amber) 163 denom = id12*id12*id32*id32 164 energy = k*(theta - theta0)**2 165 fscalar = 2.0_dp*k*(theta - theta0)/SIN(theta) 166 g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom 167 g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom 168 g3 = (b12*d12*d32 - dist*id32*d12*b32)*denom 169 CASE (do_ff_cubic) 170 denom = id12*id12*id32*id32 171 energy = k*(theta - theta0)**2*(1.0_dp + cb*(theta - theta0)) 172 fscalar = (2.0_dp*k*(theta - theta0)*(1.0_dp + cb*(theta - theta0)) + k*(theta - theta0)**2*cb)/SIN(theta) 173 g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom 174 g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom 175 g3 = (b12*d12*d32 - dist*id32*d12*b32)*denom 176 CASE (do_ff_mixed_bend_stretch) 177 178 ! 1) cubic term in theta (do_ff_cubic) 179 energy = k*(theta - theta0)**2*(1.0_dp + cb*(theta - theta0)) 180 fscalar = (2.0_dp*k*(theta - theta0)*(1.0_dp + cb*(theta - theta0)) + k*(theta - theta0)**2*cb)/SIN(theta) 181 denom = id12*id12*id32*id32 182 g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom*fscalar 183 g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom*fscalar 184 g3 = (b12*d12*d32 - dist*id32*d12*b32)*denom*fscalar 185 186 ! 2) stretch-stretch term 187 disp12 = d12 - r012 188 disp32 = d32 - r032 189 energy = energy + kss*disp12*disp32 190 g1 = g1 - kss*disp32*id12*b12 191 g2 = g2 + kss*disp32*id12*b12 192 g3 = g3 - kss*disp12*id32*b32 193 g2 = g2 + kss*disp12*id32*b32 194 195 ! 3) bend-stretch term 196 energy = energy + kbs12*disp12*(theta - theta0) + kbs32*disp32*(theta - theta0) 197 fscalar = (kbs12*disp12 + kbs32*disp32)/SIN(theta) 198 denom = id12*id12*id32*id32 199 200 ! 3a) bend part 201 g1 = g1 + (b32*d12*d32 - dist*d32*id12*b12)*denom*fscalar 202 g2 = g2 + (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom*fscalar 203 g3 = g3 + (b12*d12*d32 - dist*id32*d12*b32)*denom*fscalar 204 205 ! 3b) stretch part 206 g1 = g1 - kbs12*(theta - theta0)*id12*b12 207 g2 = g2 + kbs12*(theta - theta0)*id12*b12 208 g3 = g3 - kbs32*(theta - theta0)*id32*b32 209 g2 = g2 + kbs32*(theta - theta0)*id32*b32 210 211 ! fscalar is already included in g1, g2 and g3 212 fscalar = 1.0_dp 213 214 CASE (do_ff_harmonic, do_ff_g87) 215 denom = id12*id12*id32*id32 216 energy = f12*k*(theta - theta0)**2 217 fscalar = k*(theta - theta0)/SIN(theta) 218 g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom 219 g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom 220 g3 = (b12*d12*d32 - dist*id32*d12*b32)*denom 221 222 CASE (do_ff_mm3) 223 224 ! 1) up to sixth order in theta 225 energy = k*(theta - theta0)**2*(0.5_dp + (theta - theta0)* & 226 (-0.007_dp + (theta - theta0)*(2.8E-5_dp + (theta - theta0)* & 227 (-3.5E-7_dp + (theta - theta0)*4.5E-10_dp)))) 228 229 fscalar = k*(theta - theta0)*(1.0_dp + (theta - theta0)* & 230 (-0.021_dp + (theta - theta0)*(1.12E-4_dp + & 231 (theta - theta0)*(-1.75E-6_dp + (theta - theta0)*2.7E-9_dp))))/ & 232 SIN(theta) 233 234 denom = id12*id12*id32*id32 235 g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom*fscalar 236 g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom*fscalar 237 g3 = (b12*d12*d32 - dist*id32*d12*b32)*denom*fscalar 238 ! 2) bend-stretch term 239 disp12 = d12 - r012 240 disp32 = d32 - r032 241 energy = energy + kbs12*disp12*(theta - theta0) + kbs32*disp32*(theta - theta0) 242 fscalar = (kbs12*disp12 + kbs32*disp32)/SIN(theta) 243 denom = id12*id12*id32*id32 244 245 ! 2a) bend part 246 g1 = g1 + (b32*d12*d32 - dist*d32*id12*b12)*denom*fscalar 247 g2 = g2 + (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom*fscalar 248 g3 = g3 + (b12*d12*d32 - dist*id32*d12*b32)*denom*fscalar 249 250 ! 2b) stretch part 251 g1 = g1 - kbs12*(theta - theta0)*id12*b12 252 g2 = g2 + kbs12*(theta - theta0)*id12*b12 253 g3 = g3 - kbs32*(theta - theta0)*id32*b32 254 g2 = g2 + kbs32*(theta - theta0)*id32*b32 255 256 ! fscalar is already included in g1, g2 and g3 257 fscalar = 1.0_dp 258 CASE (do_ff_legendre) 259 ! Calculates f(x) = sum_{n=0}^{nmax} c(n) * P(n,x) 260 ! 261 ! Legendre Polynomials recursion formula: 262 ! P(n+1,x) = x*(2n+1)/(n+1) * P(n,x) - n/(n+1) * P(n-1,x) n = 1, 2,... ; P(0,x) = 1, P(1,x) = x, ... 263 ! P(n+1,x) = alpha(n,x) * P(n,x) + beta(n,x) * P(n-1,x) 264 ! with alpha(n,x) = x*(2n+1)/(n+1) 265 ! and beta(n,x) = -n/(n+1) 266 ! 267 ! Therefore 268 ! f(x) = sum_{n=0}^{nmax} c(n) * P(n,x) 269 ! can be calculated using a Clenshaw recursion 270 ! (c(n) are the coefficients for the Legendre polynomial expansion) 271 ! 272 ! f(x) = beta(1,x)*P(0,x)*y(2) + P(1,x)*y(1) + P(0,x)*c(0) 273 ! beta(1,x) = -0.5 274 ! y(k) is calculated in the following way: 275 ! y(nmax+1) = 0 276 ! y(nmax+2) = 0 277 ! y(k) = alpha(k,x)*y(k+1) + beta(k+1,x)*y(k+2) + c(k) 278 279 ! Calculates f'(x) = sum_{n=0}^{nmax} c(n) * P'(n,x) 280 ! 281 ! Recursion formula for the m-th derivatives of Legendre Polynomials 282 ! P^{m}(n+1,x) = x*(2n+1)/(n+1-m) * P^{m}(n,x) -(n+m)/(n+1-m) * P^{m}(n-1,x) n = 1, 2, ... ; m = 1, ..., n-1 283 ! For first derivative: 284 ! P'(n+1,x) = x*(2n+1)/n * P'(n,x) - (n+1)/n * P'(n-1,x) with P(0,x) = 0; P(1,x) = 1 285 ! P'(n+1,x) = alpha(n,x) * P'(n,x) + beta(n,x) * P'(n-1,x) 286 ! with alpha(n,x) = x*(2n+1)/n 287 ! and beta(n,x) = -1*(n+1)/n 288 ! 289 ! Therefore Clenshaw recursion can be used 290 ! f'(x) = beta(1,x)*P'(0,x)*y(2) + P'(1,x)*y(1) + P'(0,x)*c(0) 291 ! = beta(1,x)*0*y(2) + y(1) + 0 292 ! = y(1) 293 ! y(k) is calculated in the following way: 294 ! y(nmax+1) = 0 295 ! y(nmax+2) = 0 296 ! y(k) = alpha(k,x)*y(k+1) + beta(k+1,x)*y(k+2) + c(k) 297 ! 298 ctheta = COS(theta) 299 y1 = 0.0d0 300 y2 = 0.0d0 301 yd1 = 0.0d0 302 yd2 = 0.0d0 303 DO i = legendre%order, 2, -1 304 y0 = (2*i - 1)*ctheta*y1/i - i*y2/(i + 1) + legendre%coeffs(i) 305 y2 = y1 306 y1 = y0 307 yd0 = (2*i - 1)*ctheta*yd1/(i - 1) - (i + 1)*yd2/i + legendre%coeffs(i) 308 yd2 = yd1 309 yd1 = yd0 310 END DO 311 312 energy = -f12*y2 + ctheta*y1 + legendre%coeffs(1) 313 fscalar = -yd1 314 g1 = (b32*id32 - b12*id12*ctheta)*id12 315 g3 = (b12*id12 - b32*id32*ctheta)*id32 316 g2 = -g1 - g3 317 318 CASE DEFAULT 319 CPABORT("Unmatched bend kind") 320 END SELECT 321 322 END SUBROUTINE force_bends 323 324! ************************************************************************************************** 325!> \brief Computes the forces from the torsions 326!> \param id_type ... 327!> \param s32 ... 328!> \param is32 ... 329!> \param ism ... 330!> \param isn ... 331!> \param dist1 ... 332!> \param dist2 ... 333!> \param tm ... 334!> \param tn ... 335!> \param t12 ... 336!> \param k ... 337!> \param phi0 ... 338!> \param m ... 339!> \param gt1 ... 340!> \param gt2 ... 341!> \param gt3 ... 342!> \param gt4 ... 343!> \param energy ... 344!> \param fscalar ... 345!> \par History 346!> none 347!> \author DG 348! ************************************************************************************************** 349 SUBROUTINE force_torsions(id_type, s32, is32, ism, isn, dist1, dist2, tm, & 350 tn, t12, k, phi0, m, gt1, gt2, gt3, gt4, energy, fscalar) 351 INTEGER, INTENT(IN) :: id_type 352 REAL(KIND=dp), INTENT(IN) :: s32, is32, ism, isn, dist1, dist2 353 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tm, tn, t12 354 REAL(KIND=dp), INTENT(IN) :: k, phi0 355 INTEGER, INTENT(IN) :: m 356 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4 357 REAL(KIND=dp), INTENT(OUT) :: energy, fscalar 358 359 CHARACTER(len=*), PARAMETER :: routineN = 'force_torsions', routineP = moduleN//':'//routineN 360 361 REAL(KIND=dp) :: cosphi, phi 362 363 cosphi = DOT_PRODUCT(tm, tn)*ism*isn 364 IF (cosphi > 1.0_dp) cosphi = 1.0_dp 365 IF (cosphi < -1.0_dp) cosphi = -1.0_dp 366 phi = SIGN(ACOS(cosphi), DOT_PRODUCT(t12, tn)) 367 368 ! Select force field 369 SELECT CASE (id_type) 370 CASE (do_ff_charmm, do_ff_g87, do_ff_g96, do_ff_amber, do_ff_opls) 371 ! compute energy 372 energy = k*(1.0_dp + COS(m*phi - phi0)) 373 374 ! compute fscalar 375 fscalar = k*m*SIN(m*phi - phi0) 376 CASE DEFAULT 377 CPABORT("Unmatched torsion kind") 378 END SELECT 379 380 ! compute the gradients 381 gt1 = (s32*ism*ism)*tm 382 gt4 = -(s32*isn*isn)*tn 383 gt2 = (dist1*is32**2 - 1.0_dp)*gt1 - dist2*is32**2*gt4 384 gt3 = (dist2*is32**2 - 1.0_dp)*gt4 - dist1*is32**2*gt1 385 END SUBROUTINE force_torsions 386 387! ************************************************************************************************** 388!> \brief Computes the forces from the improper torsions 389!> \param id_type ... 390!> \param s32 ... 391!> \param is32 ... 392!> \param ism ... 393!> \param isn ... 394!> \param dist1 ... 395!> \param dist2 ... 396!> \param tm ... 397!> \param tn ... 398!> \param t12 ... 399!> \param k ... 400!> \param phi0 ... 401!> \param gt1 ... 402!> \param gt2 ... 403!> \param gt3 ... 404!> \param gt4 ... 405!> \param energy ... 406!> \param fscalar ... 407!> \par History 408!> none 409!> \author DG 410! ************************************************************************************************** 411 SUBROUTINE force_imp_torsions(id_type, s32, is32, ism, isn, dist1, dist2, tm, & 412 tn, t12, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar) 413 INTEGER, INTENT(IN) :: id_type 414 REAL(KIND=dp), INTENT(IN) :: s32, is32, ism, isn, dist1, dist2 415 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tm, tn, t12 416 REAL(KIND=dp), INTENT(IN) :: k, phi0 417 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4 418 REAL(KIND=dp), INTENT(OUT) :: energy, fscalar 419 420 CHARACTER(len=*), PARAMETER :: routineN = 'force_imp_torsions', & 421 routineP = moduleN//':'//routineN 422 REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp 423 424 REAL(KIND=dp) :: cosphi, phi 425 426! define cosphi 427 428 cosphi = DOT_PRODUCT(tm, tn)*ism*isn 429 IF (cosphi > 1.0_dp) cosphi = 1.0_dp 430 IF (cosphi < -1.0_dp) cosphi = -1.0_dp 431 phi = SIGN(ACOS(cosphi), DOT_PRODUCT(t12, tn)) 432 433 SELECT CASE (id_type) 434 CASE (do_ff_charmm) 435 ! compute energy 436 energy = k*(phi - phi0)**2 437 438 ! compute fscalar 439 fscalar = -2.0_dp*k*(phi - phi0) 440 441 CASE (do_ff_harmonic, do_ff_g87, do_ff_g96) 442 ! compute energy 443 energy = f12*k*(phi - phi0)**2 444 445 ! compute fscalar 446 fscalar = -k*(phi - phi0) 447 448 CASE DEFAULT 449 CPABORT("Unmatched improper kind") 450 END SELECT 451 452 ! compute the gradients 453 gt1 = (s32*ism*ism)*tm 454 gt4 = -(s32*isn*isn)*tn 455 gt2 = (dist1*is32**2 - 1.0_dp)*gt1 - dist2*is32**2*gt4 456 gt3 = (dist2*is32**2 - 1.0_dp)*gt4 - dist1*is32**2*gt1 457 END SUBROUTINE force_imp_torsions 458 459 ! ************************************************************************************************** 460!> \brief Computes the forces from the out of plane bends 461!> \param id_type ... 462!> \param s32 ... 463!> \param tm ... 464!> \param t41 ... 465!> \param t42 ... 466!> \param t43 ... 467!> \param k ... 468!> \param phi0 ... 469!> \param gt1 ... 470!> \param gt2 ... 471!> \param gt3 ... 472!> \param gt4 ... 473!> \param energy ... 474!> \param fscalar ... 475!> \par History 476!> none 477!> \author Louis Vanduyfhuys 478! ************************************************************************************************** 479 SUBROUTINE force_opbends(id_type, s32, tm, & 480 t41, t42, t43, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar) 481 482 INTEGER, INTENT(IN) :: id_type 483 REAL(KIND=dp), INTENT(IN) :: s32 484 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tm, t41, t42, t43 485 REAL(KIND=dp), INTENT(IN) :: k, phi0 486 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gt1, gt2, gt3, gt4 487 REAL(KIND=dp), INTENT(OUT) :: energy, fscalar 488 489 CHARACTER(len=*), PARAMETER :: routineN = 'force_opbends', routineP = moduleN//':'//routineN 490 REAL(KIND=dp), PARAMETER :: f12 = 1.0_dp/2.0_dp 491 492 REAL(KIND=dp) :: b, C, cosphi, D, E, is41, phi 493 REAL(KIND=dp), DIMENSION(3) :: db_dq1, db_dq2, db_dq3, dC_dq1, dC_dq2, & 494 dC_dq3, dD_dq1, dD_dq2, dD_dq3, & 495 dE_dq1, dE_dq2, dE_dq3 496 497!compute the energy and the gradients of cos(phi), see 498! "Efficient treatment of out-of-plane bend and improper torsion interactions in 499! MM2, MM3 and MM4 Molecular mechanicsd calculations.", Robert E. Tuzun, Donald W.Noid, 500! Bobby G Sumpter, Chemical and Analytical Sciences Division, Oak Ridge National 501! Laboratory, P.O. Box 2008, Oak Ridge, Tennesse 37831-6497 502!CAUTION: in the paper r_ij = rj - ri 503!in fist_intra_force.F t_ij = ri - rj 504!hence a minus sign needs to be added to convert r_ij vectors in t_ij vectors 505!tm is the normal of the plane 123, tm = t32 x t12 (= w from paper) 506!tn = - t41 x tm (= a from paper, for minus sign see CAUTION above) 507!Computing some necessary variables (see paper) 508 509 E = DOT_PRODUCT(-t41, tm) 510 C = DOT_PRODUCT(tm, tm) 511 D = E**2/C 512 b = DOT_PRODUCT(t41, t41) - D 513 514 !inverse norm of t41 515 is41 = 1.0_dp/SQRT(DOT_PRODUCT(t41, t41)) 516 517 cosphi = SQRT(b)*is41 518 IF (cosphi > 1.0_dp) cosphi = 1.0_dp 519 IF (cosphi < -1.0_dp) cosphi = -1.0_dp 520 phi = SIGN(ACOS(cosphi), DOT_PRODUCT(tm, -t41)) 521 522 SELECT CASE (id_type) 523 CASE (do_ff_mm2, do_ff_mm3, do_ff_mm4) 524 ! compute energy 525 energy = k*(phi - phi0)**2 526 527 ! compute fscalar 528 fscalar = 2.0_dp*k*(phi - phi0)*is41 529 530 CASE (do_ff_harmonic) 531 ! compute energy 532 energy = f12*k*(phi - phi0)**2 533 534 ! compute fscalar 535 fscalar = k*(phi - phi0)*is41 536 537 CASE DEFAULT 538 CPABORT("Unmatched opbend kind") 539 END SELECT 540 541 !Computing the necessary intermediate variables. dX_dqi is the gradient 542 !of X with respect to the coordinates of particle i. 543 544 dE_dq1(1) = (t42(2)*t43(3) - t43(2)*t42(3)) 545 dE_dq1(2) = (-t42(1)*t43(3) + t43(1)*t42(3)) 546 dE_dq1(3) = (t42(1)*t43(2) - t43(1)*t42(2)) 547 548 dE_dq2(1) = (t43(2)*t41(3) - t41(2)*t43(3)) 549 dE_dq2(2) = (-t43(1)*t41(3) + t41(1)*t43(3)) 550 dE_dq2(3) = (t43(1)*t41(2) - t41(1)*t43(2)) 551 552 dE_dq3(1) = (t41(2)*t42(3) - t42(2)*t41(3)) 553 dE_dq3(2) = (-t41(1)*t42(3) + t42(1)*t41(3)) 554 dE_dq3(3) = (t41(1)*t42(2) - t42(1)*t41(2)) 555 556 dC_dq1 = 2.0_dp*((t42 - t41)*s32**2 - (t42 - t43)*DOT_PRODUCT(t42 - t41, t42 - t43)) 557 dC_dq3 = 2.0_dp*((t42 - t43)*DOT_PRODUCT(t41 - t42, t41 - t42) & 558 - (t42 - t41)*DOT_PRODUCT(t42 - t41, t42 - t43)) 559 !C only dependent of atom 1 2 and 3, using translational invariance we find 560 dC_dq2 = -(dC_dq1 + dC_dq3) 561 562 dD_dq1 = (2.0_dp*E*dE_dq1 - D*dC_dq1)/C 563 dD_dq2 = (2.0_dp*E*dE_dq2 - D*dC_dq2)/C 564 dD_dq3 = (2.0_dp*E*dE_dq3 - D*dC_dq3)/C 565 566 db_dq1 = -2.0_dp*t41 - dD_dq1 567 db_dq2 = -dD_dq2 568 db_dq3 = -dD_dq3 569 570 !gradients of cos(phi), gt4 is calculated using translational invariance. 571 !The 1/r41 factor from the paper is absorbed in fscalar. 572 !If phi = 0 then sin(phi) = 0 and the regular formula for calculating gt 573 !won't work because of the sine function in the denominator. A further 574 !analytic simplification is needed. 575 IF (phi == 0) THEN 576 gt1 = -SIGN(1.0_dp, phi)/SQRT(C)*dE_dq1 577 gt2 = -SIGN(1.0_dp, phi)/SQRT(C)*dE_dq2 578 gt3 = -SIGN(1.0_dp, phi)/SQRT(C)*dE_dq3 579 gt4 = -(gt1 + gt2 + gt3) 580 581 ELSE 582 gt1 = (1.0_dp/(2.0_dp*SQRT(b))*db_dq1 + cosphi*t41*is41)/SIN(phi) 583 gt2 = (1.0_dp/(2.0_dp*SQRT(b))*db_dq2)/SIN(phi) 584 gt3 = (1.0_dp/(2.0_dp*SQRT(b))*db_dq3)/SIN(phi) 585 gt4 = -(gt1 + gt2 + gt3) 586 END IF 587 END SUBROUTINE force_opbends 588 589! ************************************************************************************************** 590!> \brief Computes the pressure tensor from the bonds 591!> \param f12 ... 592!> \param r12 ... 593!> \param pv_bond ... 594!> \par History 595!> none 596!> \author CJM 597! ************************************************************************************************** 598 SUBROUTINE get_pv_bond(f12, r12, pv_bond) 599 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: f12, r12 600 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pv_bond 601 602 CHARACTER(len=*), PARAMETER :: routineN = 'get_pv_bond', routineP = moduleN//':'//routineN 603 604 pv_bond(1, 1) = pv_bond(1, 1) + f12(1)*r12(1) 605 pv_bond(1, 2) = pv_bond(1, 2) + f12(1)*r12(2) 606 pv_bond(1, 3) = pv_bond(1, 3) + f12(1)*r12(3) 607 pv_bond(2, 1) = pv_bond(2, 1) + f12(2)*r12(1) 608 pv_bond(2, 2) = pv_bond(2, 2) + f12(2)*r12(2) 609 pv_bond(2, 3) = pv_bond(2, 3) + f12(2)*r12(3) 610 pv_bond(3, 1) = pv_bond(3, 1) + f12(3)*r12(1) 611 pv_bond(3, 2) = pv_bond(3, 2) + f12(3)*r12(2) 612 pv_bond(3, 3) = pv_bond(3, 3) + f12(3)*r12(3) 613 614 END SUBROUTINE get_pv_bond 615 616! ************************************************************************************************** 617!> \brief Computes the pressure tensor from the bends 618!> \param f1 ... 619!> \param f3 ... 620!> \param r12 ... 621!> \param r32 ... 622!> \param pv_bend ... 623!> \par History 624!> none 625!> \author CJM 626! ************************************************************************************************** 627 SUBROUTINE get_pv_bend(f1, f3, r12, r32, pv_bend) 628 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: f1, f3, r12, r32 629 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pv_bend 630 631 CHARACTER(len=*), PARAMETER :: routineN = 'get_pv_bend', routineP = moduleN//':'//routineN 632 633 pv_bend(1, 1) = pv_bend(1, 1) + f1(1)*r12(1) 634 pv_bend(1, 1) = pv_bend(1, 1) + f3(1)*r32(1) 635 pv_bend(1, 2) = pv_bend(1, 2) + f1(1)*r12(2) 636 pv_bend(1, 2) = pv_bend(1, 2) + f3(1)*r32(2) 637 pv_bend(1, 3) = pv_bend(1, 3) + f1(1)*r12(3) 638 pv_bend(1, 3) = pv_bend(1, 3) + f3(1)*r32(3) 639 pv_bend(2, 1) = pv_bend(2, 1) + f1(2)*r12(1) 640 pv_bend(2, 1) = pv_bend(2, 1) + f3(2)*r32(1) 641 pv_bend(2, 2) = pv_bend(2, 2) + f1(2)*r12(2) 642 pv_bend(2, 2) = pv_bend(2, 2) + f3(2)*r32(2) 643 pv_bend(2, 3) = pv_bend(2, 3) + f1(2)*r12(3) 644 pv_bend(2, 3) = pv_bend(2, 3) + f3(2)*r32(3) 645 pv_bend(3, 1) = pv_bend(3, 1) + f1(3)*r12(1) 646 pv_bend(3, 1) = pv_bend(3, 1) + f3(3)*r32(1) 647 pv_bend(3, 2) = pv_bend(3, 2) + f1(3)*r12(2) 648 pv_bend(3, 2) = pv_bend(3, 2) + f3(3)*r32(2) 649 pv_bend(3, 3) = pv_bend(3, 3) + f1(3)*r12(3) 650 pv_bend(3, 3) = pv_bend(3, 3) + f3(3)*r32(3) 651 652 END SUBROUTINE get_pv_bend 653 654! ************************************************************************************************** 655!> \brief Computes the pressure tensor from the torsions (also used for impr 656!> and opbend) 657!> \param f1 ... 658!> \param f3 ... 659!> \param f4 ... 660!> \param r12 ... 661!> \param r32 ... 662!> \param r43 ... 663!> \param pv_torsion ... 664!> \par History 665!> none 666!> \author DG 667! ************************************************************************************************** 668 SUBROUTINE get_pv_torsion(f1, f3, f4, r12, r32, r43, pv_torsion) 669 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: f1, f3, f4, r12, r32, r43 670 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pv_torsion 671 672 CHARACTER(len=*), PARAMETER :: routineN = 'get_pv_torsion', routineP = moduleN//':'//routineN 673 674 pv_torsion(1, 1) = pv_torsion(1, 1) + f1(1)*r12(1) 675 pv_torsion(1, 1) = pv_torsion(1, 1) + (f3(1) + f4(1))*r32(1) 676 pv_torsion(1, 1) = pv_torsion(1, 1) + f4(1)*r43(1) 677 pv_torsion(1, 2) = pv_torsion(1, 2) + f1(1)*r12(2) 678 pv_torsion(1, 2) = pv_torsion(1, 2) + (f3(1) + f4(1))*r32(2) 679 pv_torsion(1, 2) = pv_torsion(1, 2) + f4(1)*r43(2) 680 pv_torsion(1, 3) = pv_torsion(1, 3) + f1(1)*r12(3) 681 pv_torsion(1, 3) = pv_torsion(1, 3) + (f3(1) + f4(1))*r32(3) 682 pv_torsion(1, 3) = pv_torsion(1, 3) + f4(1)*r43(3) 683 pv_torsion(2, 1) = pv_torsion(2, 1) + f1(2)*r12(1) 684 pv_torsion(2, 1) = pv_torsion(2, 1) + (f3(2) + f4(2))*r32(1) 685 pv_torsion(2, 1) = pv_torsion(2, 1) + f4(2)*r43(1) 686 pv_torsion(2, 2) = pv_torsion(2, 2) + f1(2)*r12(2) 687 pv_torsion(2, 2) = pv_torsion(2, 2) + (f3(2) + f4(2))*r32(2) 688 pv_torsion(2, 2) = pv_torsion(2, 2) + f4(2)*r43(2) 689 pv_torsion(2, 3) = pv_torsion(2, 3) + f1(2)*r12(3) 690 pv_torsion(2, 3) = pv_torsion(2, 3) + (f3(2) + f4(2))*r32(3) 691 pv_torsion(2, 3) = pv_torsion(2, 3) + f4(2)*r43(3) 692 pv_torsion(3, 1) = pv_torsion(3, 1) + f1(3)*r12(1) 693 pv_torsion(3, 1) = pv_torsion(3, 1) + (f3(3) + f4(3))*r32(1) 694 pv_torsion(3, 1) = pv_torsion(3, 1) + f4(3)*r43(1) 695 pv_torsion(3, 2) = pv_torsion(3, 2) + f1(3)*r12(2) 696 pv_torsion(3, 2) = pv_torsion(3, 2) + (f3(3) + f4(3))*r32(2) 697 pv_torsion(3, 2) = pv_torsion(3, 2) + f4(3)*r43(2) 698 pv_torsion(3, 3) = pv_torsion(3, 3) + f1(3)*r12(3) 699 pv_torsion(3, 3) = pv_torsion(3, 3) + (f3(3) + f4(3))*r32(3) 700 pv_torsion(3, 3) = pv_torsion(3, 3) + f4(3)*r43(3) 701 702 END SUBROUTINE get_pv_torsion 703 704END MODULE mol_force 705 706