(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE 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
29CONTAINS
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
78 CASE (do_ff_charmm, do_ff_amber)
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
88 CASE (do_ff_harmonic, do_ff_g87)
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
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 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)
368 CASE (do_ff_charmm, do_ff_g87, do_ff_g96, do_ff_amber, do_ff_opls)
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
437 CASE (do_ff_harmonic, do_ff_g87, do_ff_g96)
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
693END MODULE mol_force
694
force_field_kind_types
Define all structure types related to force field kinds.
Definition force_field_kind_types.F:14
force_field_kind_types::do_ff_legendre
integer, parameter, public do_ff_legendre
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_mm4
integer, parameter, public do_ff_mm4
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_charmm
integer, parameter, public do_ff_charmm
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_mm3
integer, parameter, public do_ff_mm3
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_g87
integer, parameter, public do_ff_g87
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_g96
integer, parameter, public do_ff_g96
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_morse
integer, parameter, public do_ff_morse
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_mm2
integer, parameter, public do_ff_mm2
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_harmonic
integer, parameter, public do_ff_harmonic
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_amber
integer, parameter, public do_ff_amber
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_mixed_bend_stretch
integer, parameter, public do_ff_mixed_bend_stretch
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_cubic
integer, parameter, public do_ff_cubic
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_quartic
integer, parameter, public do_ff_quartic
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_fues
integer, parameter, public do_ff_fues
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_opls
integer, parameter, public do_ff_opls
Definition force_field_kind_types.F:25
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mol_force
Definition mol_force.F:13
mol_force::force_opbends
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
mol_force::force_imp_torsions
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
mol_force::get_pv_bend
subroutine, public get_pv_bend(f1, f3, r12, r32, pv_bend)
Computes the pressure tensor from the bends.
Definition mol_force.F:621
mol_force::get_pv_bond
subroutine, public get_pv_bond(f12, r12, pv_bond)
Computes the pressure tensor from the bonds.
Definition mol_force.F:594
mol_force::force_bonds
subroutine, public force_bonds(id_type, rij, r0, k, cs, energy, fscalar)
Computes the forces from the bonds.
Definition mol_force.F:43
mol_force::get_pv_torsion
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
mol_force::force_torsions
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
mol_force::force_bends
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
force_field_kind_types::legendre_data_type
Definition force_field_kind_types.F:45