LCOV - code coverage report
Current view: top level - src - mol_force.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.7 % 263 257
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 8 8

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

Generated by: LCOV version 2.0-1