(git:0de0cc2)
mol_force.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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: &
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'
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  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  SELECT CASE (id_type)
55  CASE (do_ff_quartic)
56  dij = sqrt(dot_product(rij, rij))
57  disp = dij - r0
58  energy = (f12*k(1) + (f13*k(2) + f14*k(3)*disp)*disp)*disp*disp
59  fscalar = ((k(1) + (k(2) + k(3)*disp)*disp)*disp)/dij
60  CASE (do_ff_morse)
61  dij = sqrt(dot_product(rij, rij))
62  disp = dij - r0
63  energy = k(1)*((1 - exp(-k(2)*disp))**2 - 1)
64  fscalar = 2*k(1)*k(2)*exp(-k(2)*disp)*(1 - exp(-k(2)*disp))/dij
65  CASE (do_ff_cubic)
66  dij = sqrt(dot_product(rij, rij))
67  disp = dij - r0
68  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  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  dij = dot_product(rij, rij)
75  disp = dij - r0*r0
76  energy = f14*k(1)*disp*disp
77  fscalar = k(1)*disp
79  dij = sqrt(dot_product(rij, rij))
80  disp = dij - r0
81  IF (abs(disp) < epsilon(1.0_dp)) THEN
82  energy = 0.0_dp
83  fscalar = 0.0_dp
84  ELSE
85  energy = k(1)*disp*disp
86  fscalar = 2.0_dp*k(1)*disp/dij
87  END IF
89  dij = sqrt(dot_product(rij, rij))
90  disp = dij - r0
91  IF (abs(disp) < epsilon(1.0_dp)) THEN
92  energy = 0.0_dp
93  fscalar = 0.0_dp
94  ELSE
95  energy = f12*k(1)*disp*disp
96  fscalar = k(1)*disp/dij
97  END IF
98  CASE (do_ff_fues)
99  dij = sqrt(dot_product(rij, rij))
100  disp = r0/dij
101  energy = f12*k(1)*r0*r0*(1.0_dp + disp*(disp - 2.0_dp))
102  fscalar = k(1)*r0*disp*disp*(1.0_dp - disp)/dij
103  CASE DEFAULT
104  cpabort("Unmatched bond kind")
105  END SELECT
106 
107  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  SUBROUTINE force_bends(id_type, b12, b32, d12, d32, id12, id32, dist, &
139  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  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
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  REAL(kind=dp) :: cosphi, phi
360 
361  cosphi = dot_product(tm, tn)*ism*isn
362  IF (cosphi > 1.0_dp) cosphi = 1.0_dp
363  IF (cosphi < -1.0_dp) cosphi = -1.0_dp
364  phi = sign(acos(cosphi), dot_product(t12, tn))
365 
366  ! Select force field
367  SELECT CASE (id_type)
369  ! compute energy
370  energy = k*(1.0_dp + cos(m*phi - phi0))
371 
372  ! compute fscalar
373  fscalar = k*m*sin(m*phi - phi0)
374  CASE DEFAULT
375  cpabort("Unmatched torsion kind")
376  END SELECT
377 
378  ! compute the gradients
379  gt1 = (s32*ism*ism)*tm
380  gt4 = -(s32*isn*isn)*tn
381  gt2 = (dist1*is32**2 - 1.0_dp)*gt1 - dist2*is32**2*gt4
382  gt3 = (dist2*is32**2 - 1.0_dp)*gt4 - dist1*is32**2*gt1
383  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  SUBROUTINE force_imp_torsions(id_type, s32, is32, ism, isn, dist1, dist2, tm, &
410  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  cosphi = dot_product(tm, tn)*ism*isn
425  IF (cosphi > 1.0_dp) cosphi = 1.0_dp
426  IF (cosphi < -1.0_dp) cosphi = -1.0_dp
427  phi = sign(acos(cosphi), dot_product(t12, tn))
428 
429  SELECT CASE (id_type)
430  CASE (do_ff_charmm)
431  ! compute energy
432  energy = k*(phi - phi0)**2
433 
434  ! compute fscalar
435  fscalar = -2.0_dp*k*(phi - phi0)
436 
438  ! compute energy
439  energy = f12*k*(phi - phi0)**2
440 
441  ! compute fscalar
442  fscalar = -k*(phi - phi0)
443 
444  CASE DEFAULT
445  cpabort("Unmatched improper kind")
446  END SELECT
447 
448  ! compute the gradients
449  gt1 = (s32*ism*ism)*tm
450  gt4 = -(s32*isn*isn)*tn
451  gt2 = (dist1*is32**2 - 1.0_dp)*gt1 - dist2*is32**2*gt4
452  gt3 = (dist2*is32**2 - 1.0_dp)*gt4 - dist1*is32**2*gt1
453  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  SUBROUTINE force_opbends(id_type, s32, tm, &
476  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  e = dot_product(-t41, tm)
505  c = dot_product(tm, tm)
506  d = e**2/c
507  b = dot_product(t41, t41) - d
508 
509  !inverse norm of t41
510  is41 = 1.0_dp/sqrt(dot_product(t41, t41))
511 
512  cosphi = sqrt(b)*is41
513  IF (cosphi > 1.0_dp) cosphi = 1.0_dp
514  IF (cosphi < -1.0_dp) cosphi = -1.0_dp
515  phi = sign(acos(cosphi), dot_product(tm, -t41))
516 
517  SELECT CASE (id_type)
518  CASE (do_ff_mm2, do_ff_mm3, do_ff_mm4)
519  ! compute energy
520  energy = k*(phi - phi0)**2
521 
522  ! compute fscalar
523  fscalar = 2.0_dp*k*(phi - phi0)*is41
524 
525  CASE (do_ff_harmonic)
526  ! compute energy
527  energy = f12*k*(phi - phi0)**2
528 
529  ! compute fscalar
530  fscalar = k*(phi - phi0)*is41
531 
532  CASE DEFAULT
533  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  de_dq1(1) = (t42(2)*t43(3) - t43(2)*t42(3))
540  de_dq1(2) = (-t42(1)*t43(3) + t43(1)*t42(3))
541  de_dq1(3) = (t42(1)*t43(2) - t43(1)*t42(2))
542 
543  de_dq2(1) = (t43(2)*t41(3) - t41(2)*t43(3))
544  de_dq2(2) = (-t43(1)*t41(3) + t41(1)*t43(3))
545  de_dq2(3) = (t43(1)*t41(2) - t41(1)*t43(2))
546 
547  de_dq3(1) = (t41(2)*t42(3) - t42(2)*t41(3))
548  de_dq3(2) = (-t41(1)*t42(3) + t42(1)*t41(3))
549  de_dq3(3) = (t41(1)*t42(2) - t42(1)*t41(2))
550 
551  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  - (t42 - t41)*dot_product(t42 - t41, t42 - t43))
554  !C only dependent of atom 1 2 and 3, using translational invariance we find
555  dc_dq2 = -(dc_dq1 + dc_dq3)
556 
557  dd_dq1 = (2.0_dp*e*de_dq1 - d*dc_dq1)/c
558  dd_dq2 = (2.0_dp*e*de_dq2 - d*dc_dq2)/c
559  dd_dq3 = (2.0_dp*e*de_dq3 - d*dc_dq3)/c
560 
561  db_dq1 = -2.0_dp*t41 - dd_dq1
562  db_dq2 = -dd_dq2
563  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  IF (phi == 0) THEN
571  gt1 = -sign(1.0_dp, phi)/sqrt(c)*de_dq1
572  gt2 = -sign(1.0_dp, phi)/sqrt(c)*de_dq2
573  gt3 = -sign(1.0_dp, phi)/sqrt(c)*de_dq3
574  gt4 = -(gt1 + gt2 + gt3)
575 
576  ELSE
577  gt1 = (1.0_dp/(2.0_dp*sqrt(b))*db_dq1 + cosphi*t41*is41)/sin(phi)
578  gt2 = (1.0_dp/(2.0_dp*sqrt(b))*db_dq2)/sin(phi)
579  gt3 = (1.0_dp/(2.0_dp*sqrt(b))*db_dq3)/sin(phi)
580  gt4 = -(gt1 + gt2 + gt3)
581  END IF
582  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  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  pv_bond(1, 1) = pv_bond(1, 1) + f12(1)*r12(1)
598  pv_bond(1, 2) = pv_bond(1, 2) + f12(1)*r12(2)
599  pv_bond(1, 3) = pv_bond(1, 3) + f12(1)*r12(3)
600  pv_bond(2, 1) = pv_bond(2, 1) + f12(2)*r12(1)
601  pv_bond(2, 2) = pv_bond(2, 2) + f12(2)*r12(2)
602  pv_bond(2, 3) = pv_bond(2, 3) + f12(2)*r12(3)
603  pv_bond(3, 1) = pv_bond(3, 1) + f12(3)*r12(1)
604  pv_bond(3, 2) = pv_bond(3, 2) + f12(3)*r12(2)
605  pv_bond(3, 3) = pv_bond(3, 3) + f12(3)*r12(3)
606 
607  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  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  pv_bend(1, 1) = pv_bend(1, 1) + f1(1)*r12(1)
625  pv_bend(1, 1) = pv_bend(1, 1) + f3(1)*r32(1)
626  pv_bend(1, 2) = pv_bend(1, 2) + f1(1)*r12(2)
627  pv_bend(1, 2) = pv_bend(1, 2) + f3(1)*r32(2)
628  pv_bend(1, 3) = pv_bend(1, 3) + f1(1)*r12(3)
629  pv_bend(1, 3) = pv_bend(1, 3) + f3(1)*r32(3)
630  pv_bend(2, 1) = pv_bend(2, 1) + f1(2)*r12(1)
631  pv_bend(2, 1) = pv_bend(2, 1) + f3(2)*r32(1)
632  pv_bend(2, 2) = pv_bend(2, 2) + f1(2)*r12(2)
633  pv_bend(2, 2) = pv_bend(2, 2) + f3(2)*r32(2)
634  pv_bend(2, 3) = pv_bend(2, 3) + f1(2)*r12(3)
635  pv_bend(2, 3) = pv_bend(2, 3) + f3(2)*r32(3)
636  pv_bend(3, 1) = pv_bend(3, 1) + f1(3)*r12(1)
637  pv_bend(3, 1) = pv_bend(3, 1) + f3(3)*r32(1)
638  pv_bend(3, 2) = pv_bend(3, 2) + f1(3)*r12(2)
639  pv_bend(3, 2) = pv_bend(3, 2) + f3(3)*r32(2)
640  pv_bend(3, 3) = pv_bend(3, 3) + f1(3)*r12(3)
641  pv_bend(3, 3) = pv_bend(3, 3) + f3(3)*r32(3)
642 
643  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  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  pv_torsion(1, 1) = pv_torsion(1, 1) + f1(1)*r12(1)
664  pv_torsion(1, 1) = pv_torsion(1, 1) + (f3(1) + f4(1))*r32(1)
665  pv_torsion(1, 1) = pv_torsion(1, 1) + f4(1)*r43(1)
666  pv_torsion(1, 2) = pv_torsion(1, 2) + f1(1)*r12(2)
667  pv_torsion(1, 2) = pv_torsion(1, 2) + (f3(1) + f4(1))*r32(2)
668  pv_torsion(1, 2) = pv_torsion(1, 2) + f4(1)*r43(2)
669  pv_torsion(1, 3) = pv_torsion(1, 3) + f1(1)*r12(3)
670  pv_torsion(1, 3) = pv_torsion(1, 3) + (f3(1) + f4(1))*r32(3)
671  pv_torsion(1, 3) = pv_torsion(1, 3) + f4(1)*r43(3)
672  pv_torsion(2, 1) = pv_torsion(2, 1) + f1(2)*r12(1)
673  pv_torsion(2, 1) = pv_torsion(2, 1) + (f3(2) + f4(2))*r32(1)
674  pv_torsion(2, 1) = pv_torsion(2, 1) + f4(2)*r43(1)
675  pv_torsion(2, 2) = pv_torsion(2, 2) + f1(2)*r12(2)
676  pv_torsion(2, 2) = pv_torsion(2, 2) + (f3(2) + f4(2))*r32(2)
677  pv_torsion(2, 2) = pv_torsion(2, 2) + f4(2)*r43(2)
678  pv_torsion(2, 3) = pv_torsion(2, 3) + f1(2)*r12(3)
679  pv_torsion(2, 3) = pv_torsion(2, 3) + (f3(2) + f4(2))*r32(3)
680  pv_torsion(2, 3) = pv_torsion(2, 3) + f4(2)*r43(3)
681  pv_torsion(3, 1) = pv_torsion(3, 1) + f1(3)*r12(1)
682  pv_torsion(3, 1) = pv_torsion(3, 1) + (f3(3) + f4(3))*r32(1)
683  pv_torsion(3, 1) = pv_torsion(3, 1) + f4(3)*r43(1)
684  pv_torsion(3, 2) = pv_torsion(3, 2) + f1(3)*r12(2)
685  pv_torsion(3, 2) = pv_torsion(3, 2) + (f3(3) + f4(3))*r32(2)
686  pv_torsion(3, 2) = pv_torsion(3, 2) + f4(3)*r43(2)
687  pv_torsion(3, 3) = pv_torsion(3, 3) + f1(3)*r12(3)
688  pv_torsion(3, 3) = pv_torsion(3, 3) + (f3(3) + f4(3))*r32(3)
689  pv_torsion(3, 3) = pv_torsion(3, 3) + f4(3)*r43(3)
690 
691  END SUBROUTINE get_pv_torsion
692 
693 END MODULE mol_force
694 
Define all structure types related to force field kinds.
integer, parameter, public do_ff_legendre
integer, parameter, public do_ff_mm4
integer, parameter, public do_ff_charmm
integer, parameter, public do_ff_mm3
integer, parameter, public do_ff_g87
integer, parameter, public do_ff_g96
integer, parameter, public do_ff_morse
integer, parameter, public do_ff_mm2
integer, parameter, public do_ff_harmonic
integer, parameter, public do_ff_amber
integer, parameter, public do_ff_mixed_bend_stretch
integer, parameter, public do_ff_cubic
integer, parameter, public do_ff_quartic
integer, parameter, public do_ff_fues
integer, parameter, public do_ff_opls
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
subroutine, public force_opbends(id_type, s32, tm, t41, t42, t43, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar)
Computes the forces from the out of plane bends.
Definition: mol_force.F:477
subroutine, public force_imp_torsions(id_type, s32, is32, ism, isn, dist1, dist2, tm, tn, t12, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar)
Computes the forces from the improper torsions.
Definition: mol_force.F:411
subroutine, public get_pv_bend(f1, f3, r12, r32, pv_bend)
Computes the pressure tensor from the bends.
Definition: mol_force.F:621
subroutine, public get_pv_bond(f12, r12, pv_bond)
Computes the pressure tensor from the bonds.
Definition: mol_force.F:594
subroutine, public force_bonds(id_type, rij, r0, k, cs, energy, fscalar)
Computes the forces from the bonds.
Definition: mol_force.F:43
subroutine, public get_pv_torsion(f1, f3, f4, r12, r32, r43, pv_torsion)
Computes the pressure tensor from the torsions (also used for impr and opbend)
Definition: mol_force.F:660
subroutine, public force_torsions(id_type, s32, is32, ism, isn, dist1, dist2, tm, tn, t12, k, phi0, m, gt1, gt2, gt3, gt4, energy, fscalar)
Computes the forces from the torsions.
Definition: mol_force.F:351
subroutine, public force_bends(id_type, b12, b32, d12, d32, id12, id32, dist, theta, theta0, k, cb, r012, r032, kbs12, kbs32, kss, legendre, g1, g2, g3, energy, fscalar)
Computes the forces from the bends.
Definition: mol_force.F:140