Line data Source code
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: &
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 2881621 : 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 2881621 : 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 2016 : dij = SQRT(DOT_PRODUCT(rij, rij))
62 504 : disp = dij - r0
63 504 : energy = k(1)*((1 - EXP(-k(2)*disp))**2 - 1)
64 504 : 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 10911196 : dij = SQRT(DOT_PRODUCT(rij, rij))
80 2727799 : disp = dij - r0
81 2727799 : IF (ABS(disp) < EPSILON(1.0_dp)) THEN
82 1 : energy = 0.0_dp
83 1 : fscalar = 0.0_dp
84 : ELSE
85 2727798 : energy = k(1)*disp*disp
86 2727798 : 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 2881621 : CPABORT("Unmatched bond kind")
105 : END SELECT
106 :
107 2881621 : 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 1351317 : SUBROUTINE force_bends(id_type, b12, b32, d12, d32, id12, id32, dist, &
139 2702634 : 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 1355357 : 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 1235515 : denom = id12*id12*id32*id32
164 1235515 : energy = k*(theta - theta0)**2
165 1235515 : fscalar = 2.0_dp*k*(theta - theta0)/SIN(theta)
166 4942060 : g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom
167 4942060 : g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom
168 4942060 : 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 111737 : denom = id12*id12*id32*id32
216 111737 : energy = f12*k*(theta - theta0)**2
217 111737 : fscalar = k*(theta - theta0)/SIN(theta)
218 446948 : g1 = (b32*d12*d32 - dist*d32*id12*b12)*denom
219 446948 : g2 = (-(b12 + b32)*d12*d32 + dist*(d32*id12*b12 + id32*d12*b32))*denom
220 446948 : 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 1351317 : CPABORT("Unmatched bend kind")
320 : END SELECT
321 :
322 1351317 : 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 541716 : 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 541716 : pv_bond(1, 1) = pv_bond(1, 1) + f12(1)*r12(1)
598 541716 : pv_bond(1, 2) = pv_bond(1, 2) + f12(1)*r12(2)
599 541716 : pv_bond(1, 3) = pv_bond(1, 3) + f12(1)*r12(3)
600 541716 : pv_bond(2, 1) = pv_bond(2, 1) + f12(2)*r12(1)
601 541716 : pv_bond(2, 2) = pv_bond(2, 2) + f12(2)*r12(2)
602 541716 : pv_bond(2, 3) = pv_bond(2, 3) + f12(2)*r12(3)
603 541716 : pv_bond(3, 1) = pv_bond(3, 1) + f12(3)*r12(1)
604 541716 : pv_bond(3, 2) = pv_bond(3, 2) + f12(3)*r12(2)
605 541716 : pv_bond(3, 3) = pv_bond(3, 3) + f12(3)*r12(3)
606 :
607 541716 : 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 110450 : 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 110450 : pv_bend(1, 1) = pv_bend(1, 1) + f1(1)*r12(1)
625 110450 : pv_bend(1, 1) = pv_bend(1, 1) + f3(1)*r32(1)
626 110450 : pv_bend(1, 2) = pv_bend(1, 2) + f1(1)*r12(2)
627 110450 : pv_bend(1, 2) = pv_bend(1, 2) + f3(1)*r32(2)
628 110450 : pv_bend(1, 3) = pv_bend(1, 3) + f1(1)*r12(3)
629 110450 : pv_bend(1, 3) = pv_bend(1, 3) + f3(1)*r32(3)
630 110450 : pv_bend(2, 1) = pv_bend(2, 1) + f1(2)*r12(1)
631 110450 : pv_bend(2, 1) = pv_bend(2, 1) + f3(2)*r32(1)
632 110450 : pv_bend(2, 2) = pv_bend(2, 2) + f1(2)*r12(2)
633 110450 : pv_bend(2, 2) = pv_bend(2, 2) + f3(2)*r32(2)
634 110450 : pv_bend(2, 3) = pv_bend(2, 3) + f1(2)*r12(3)
635 110450 : pv_bend(2, 3) = pv_bend(2, 3) + f3(2)*r32(3)
636 110450 : pv_bend(3, 1) = pv_bend(3, 1) + f1(3)*r12(1)
637 110450 : pv_bend(3, 1) = pv_bend(3, 1) + f3(3)*r32(1)
638 110450 : pv_bend(3, 2) = pv_bend(3, 2) + f1(3)*r12(2)
639 110450 : pv_bend(3, 2) = pv_bend(3, 2) + f3(3)*r32(2)
640 110450 : pv_bend(3, 3) = pv_bend(3, 3) + f1(3)*r12(3)
641 110450 : pv_bend(3, 3) = pv_bend(3, 3) + f3(3)*r32(3)
642 :
643 110450 : 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 :
|
