20 #include "./base/base_uses.f90"
25 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mol_force'
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
48 REAL(kind=
dp),
PARAMETER :: f12 = 1.0_dp/2.0_dp, &
49 f13 = 1.0_dp/3.0_dp, &
52 REAL(kind=
dp) :: dij, disp
56 dij = sqrt(dot_product(rij, rij))
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
61 dij = sqrt(dot_product(rij, rij))
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
66 dij = sqrt(dot_product(rij, rij))
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
74 dij = dot_product(rij, rij)
76 energy = f14*k(1)*disp*disp
79 dij = sqrt(dot_product(rij, rij))
81 IF (abs(disp) < epsilon(1.0_dp))
THEN
85 energy = k(1)*disp*disp
86 fscalar = 2.0_dp*k(1)*disp/dij
89 dij = sqrt(dot_product(rij, rij))
91 IF (abs(disp) < epsilon(1.0_dp))
THEN
95 energy = f12*k(1)*disp*disp
96 fscalar = k(1)*disp/dij
99 dij = sqrt(dot_product(rij, rij))
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
104 cpabort(
"Unmatched bond kind")
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, &
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
149 REAL(kind=
dp),
PARAMETER :: f12 = 1.0_dp/2.0_dp
152 REAL(kind=
dp) :: ctheta, denom, disp12, disp32, y0, y1, &
155 SELECT CASE (id_type)
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
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
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
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
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
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
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
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
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
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))))
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))))/ &
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
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
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
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
303 DO i = legendre%order, 2, -1
304 y0 = (2*i - 1)*ctheta*y1/i - i*y2/(i + 1) + legendre%coeffs(i)
307 yd0 = (2*i - 1)*ctheta*yd1/(i - 1) - (i + 1)*yd2/i + legendre%coeffs(i)
312 energy = -f12*y2 + ctheta*y1 + legendre%coeffs(1)
314 g1 = (b32*id32 - b12*id12*ctheta)*id12
315 g3 = (b12*id12 - b32*id32*ctheta)*id32
319 cpabort(
"Unmatched bend kind")
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
359 REAL(kind=
dp) :: cosphi, phi
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))
367 SELECT CASE (id_type)
370 energy = k*(1.0_dp + cos(m*phi - phi0))
373 fscalar = k*m*sin(m*phi - phi0)
375 cpabort(
"Unmatched torsion kind")
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
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
418 REAL(kind=
dp),
PARAMETER :: f12 = 1.0_dp/2.0_dp
420 REAL(kind=
dp) :: cosphi, phi
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))
429 SELECT CASE (id_type)
432 energy = k*(phi - phi0)**2
435 fscalar = -2.0_dp*k*(phi - phi0)
439 energy = f12*k*(phi - phi0)**2
442 fscalar = -k*(phi - phi0)
445 cpabort(
"Unmatched improper kind")
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
476 t41, t42, t43, k, phi0, gt1, gt2, gt3, gt4, energy, fscalar)
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
485 REAL(kind=
dp),
PARAMETER :: f12 = 1.0_dp/2.0_dp
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
504 e = dot_product(-t41, tm)
505 c = dot_product(tm, tm)
507 b = dot_product(t41, t41) - d
510 is41 = 1.0_dp/sqrt(dot_product(t41, t41))
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))
517 SELECT CASE (id_type)
520 energy = k*(phi - phi0)**2
523 fscalar = 2.0_dp*k*(phi - phi0)*is41
527 energy = f12*k*(phi - phi0)**2
530 fscalar = k*(phi - phi0)*is41
533 cpabort(
"Unmatched opbend kind")
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))
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))
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))
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))
555 dc_dq2 = -(dc_dq1 + dc_dq3)
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
561 db_dq1 = -2.0_dp*t41 - dd_dq1
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)
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)
594 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: f12, r12
595 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: pv_bond
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)
621 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: f1, f3, r12, r32
622 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: pv_bend
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)
660 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: f1, f3, f4, r12, r32, r43
661 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: pv_torsion
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)
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.
integer, parameter, public dp
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.
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.
subroutine, public get_pv_bend(f1, f3, r12, r32, pv_bend)
Computes the pressure tensor from the bends.
subroutine, public get_pv_bond(f12, r12, pv_bond)
Computes the pressure tensor from the bonds.
subroutine, public force_bonds(id_type, rij, r0, k, cs, energy, fscalar)
Computes the forces from the bonds.
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)
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.
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.