(git:374b731)
Loading...
Searching...
No Matches
fist_intra_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!> Variable names changed (DG) 05-Dec-2000
12!> \author CJM
13! **************************************************************************************************
14MODULE fist_intra_force
15
16 USE atprop_types, ONLY: atprop_type
17 USE cell_types, ONLY: cell_type,&
18 pbc
19 USE cp_log_handling, ONLY: cp_get_default_logger,&
20 cp_logger_type
21 USE distribution_1d_types, ONLY: distribution_1d_type
22 USE kinds, ONLY: dp
23 USE mol_force, ONLY: force_bends,&
24 force_bonds,&
25 force_imp_torsions,&
26 force_opbends,&
27 force_torsions,&
28 get_pv_bend,&
29 get_pv_bond,&
30 get_pv_torsion
31 USE molecule_kind_types, ONLY: &
32 bend_type, bond_type, get_molecule_kind, impr_type, molecule_kind_type, opbend_type, &
33 shell_type, torsion_type, ub_type
34 USE molecule_types, ONLY: get_molecule,&
35 molecule_type
36 USE particle_types, ONLY: particle_type
37#include "./base/base_uses.f90"
38
39 IMPLICIT NONE
40
41 PRIVATE
42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_intra_force'
43 PUBLIC :: force_intra_control
44
45CONTAINS
46
47! **************************************************************************************************
48!> \brief Computes the the intramolecular energies, forces, and pressure tensors
49!> \param molecule_set ...
50!> \param molecule_kind_set ...
51!> \param local_molecules ...
52!> \param particle_set ...
53!> \param shell_particle_set ...
54!> \param core_particle_set ...
55!> \param pot_bond ...
56!> \param pot_bend ...
57!> \param pot_urey_bradley ...
58!> \param pot_torsion ...
59!> \param pot_imp_torsion ...
60!> \param pot_opbend ...
61!> \param pot_shell ...
62!> \param pv_bond ...
63!> \param pv_bend ...
64!> \param pv_urey_bradley ...
65!> \param pv_torsion ...
66!> \param pv_imp_torsion ...
67!> \param pv_opbend ...
68!> \param f_bond ...
69!> \param f_bend ...
70!> \param f_torsion ...
71!> \param f_ub ...
72!> \param f_imptor ...
73!> \param f_opbend ...
74!> \param cell ...
75!> \param use_virial ...
76!> \param atprop_env ...
77!> \par History
78!> none
79!> \author CJM
80! **************************************************************************************************
81 SUBROUTINE force_intra_control(molecule_set, molecule_kind_set, &
82 local_molecules, particle_set, shell_particle_set, core_particle_set, &
83 pot_bond, pot_bend, pot_urey_bradley, pot_torsion, pot_imp_torsion, &
84 pot_opbend, pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, &
85 pv_imp_torsion, pv_opbend, f_bond, f_bend, f_torsion, f_ub, &
86 f_imptor, f_opbend, cell, use_virial, atprop_env)
87
88 TYPE(molecule_type), POINTER :: molecule_set(:)
89 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
90 TYPE(distribution_1d_type), POINTER :: local_molecules
91 TYPE(particle_type), POINTER :: particle_set(:), shell_particle_set(:), &
92 core_particle_set(:)
93 REAL(kind=dp), INTENT(INOUT) :: pot_bond, pot_bend, pot_urey_bradley, &
94 pot_torsion, pot_imp_torsion, &
95 pot_opbend, pot_shell
96 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: pv_bond, pv_bend, pv_urey_bradley, &
97 pv_torsion, pv_imp_torsion, pv_opbend
98 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
99 OPTIONAL :: f_bond, f_bend, f_torsion, f_ub, &
100 f_imptor, f_opbend
101 TYPE(cell_type), POINTER :: cell
102 LOGICAL, INTENT(IN) :: use_virial
103 TYPE(atprop_type), POINTER :: atprop_env
104
105 CHARACTER(len=*), PARAMETER :: routinen = 'force_intra_control'
106
107 INTEGER :: first_atom, handle, i, ibend, ibond, ikind, imol, imul, index_a, index_b, &
108 index_c, index_d, iopbend, ishell, itorsion, nbends, nbonds, nimptors, nkind, &
109 nmol_per_kind, nopbends, nshell, ntorsions, nub
110 LOGICAL :: atener, atstress
111 REAL(kind=dp) :: d12, d32, dist, dist1, dist2, energy, &
112 fscalar, id12, id32, is32, ism, isn, &
113 k2_spring, k4_spring, r2, s32, sm, sn, &
114 theta
115 REAL(kind=dp), DIMENSION(3) :: b12, b32, g1, g2, g3, gt1, gt2, gt3, &
116 gt4, k1, k2, k3, k4, rij, t12, t32, &
117 t34, t41, t42, t43, tm, tn
118 REAL(kind=dp), DIMENSION(:), POINTER :: ener_a
119 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pv_a
120 TYPE(bend_type), POINTER :: bend_list(:)
121 TYPE(bond_type), POINTER :: bond_list(:)
122 TYPE(cp_logger_type), POINTER :: logger
123 TYPE(impr_type), POINTER :: impr_list(:)
124 TYPE(molecule_kind_type), POINTER :: molecule_kind
125 TYPE(molecule_type), POINTER :: molecule
126 TYPE(opbend_type), POINTER :: opbend_list(:)
127 TYPE(shell_type), POINTER :: shell_list(:)
128 TYPE(torsion_type), POINTER :: torsion_list(:)
129 TYPE(ub_type), POINTER :: ub_list(:)
130
131 CALL timeset(routinen, handle)
132 NULLIFY (logger)
133 logger => cp_get_default_logger()
134
135 IF (PRESENT(f_bond)) f_bond = 0.0_dp
136 IF (PRESENT(f_bend)) f_bend = 0.0_dp
137 IF (PRESENT(f_torsion)) f_torsion = 0.0_dp
138 IF (PRESENT(f_imptor)) f_imptor = 0.0_dp
139 IF (PRESENT(f_ub)) f_ub = 0.0_dp
140
141 pot_bond = 0.0_dp
142 pot_bend = 0.0_dp
143 pot_urey_bradley = 0.0_dp
144 pot_torsion = 0.0_dp
145 pot_imp_torsion = 0.0_dp
146 pot_opbend = 0.0_dp
147 pot_shell = 0.0_dp
148
149 atener = atprop_env%energy
150 IF (atener) ener_a => atprop_env%atener
151 atstress = atprop_env%stress
152 IF (atstress) pv_a => atprop_env%atstress
153
154 nkind = SIZE(molecule_kind_set)
155 mol: DO ikind = 1, nkind
156 nmol_per_kind = local_molecules%n_el(ikind)
157
158 DO imol = 1, nmol_per_kind
159 i = local_molecules%list(ikind)%array(imol)
160 molecule => molecule_set(i)
161 molecule_kind => molecule%molecule_kind
162 CALL get_molecule_kind(molecule_kind, nbend=nbends, nbond=nbonds, &
163 nimpr=nimptors, nub=nub, ntorsion=ntorsions, &
164 nopbend=nopbends, nshell=nshell, &
165 bond_list=bond_list, ub_list=ub_list, &
166 bend_list=bend_list, torsion_list=torsion_list, &
167 impr_list=impr_list, opbend_list=opbend_list, &
168 shell_list=shell_list)
169
170 CALL get_molecule(molecule, first_atom=first_atom)
171
172 bond: DO ibond = 1, nbonds
173 index_a = bond_list(ibond)%a + first_atom - 1
174 index_b = bond_list(ibond)%b + first_atom - 1
175 rij = particle_set(index_a)%r - particle_set(index_b)%r
176 rij = pbc(rij, cell)
177 CALL force_bonds(bond_list(ibond)%bond_kind%id_type, rij, &
178 bond_list(ibond)%bond_kind%r0, &
179 bond_list(ibond)%bond_kind%k, &
180 bond_list(ibond)%bond_kind%cs, &
181 energy, fscalar)
182 pot_bond = pot_bond + energy
183 IF (atener) THEN
184 ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
185 ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
186 END IF
187
188 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
189 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
190 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
191 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
192 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
193 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
194
195 ! computing the pressure tensor
196 k2 = -rij*fscalar
197 IF (use_virial) CALL get_pv_bond(k2, rij, pv_bond)
198 IF (atstress) THEN
199 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_a))
200 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_b))
201 END IF
202
203 ! the contribution from the bonds. ONLY FOR DEBUG
204 IF (PRESENT(f_bond)) THEN
205 f_bond(1, index_a) = f_bond(1, index_a) - rij(1)*fscalar
206 f_bond(2, index_a) = f_bond(2, index_a) - rij(2)*fscalar
207 f_bond(3, index_a) = f_bond(3, index_a) - rij(3)*fscalar
208 f_bond(1, index_b) = f_bond(1, index_b) + rij(1)*fscalar
209 f_bond(2, index_b) = f_bond(2, index_b) + rij(2)*fscalar
210 f_bond(3, index_b) = f_bond(3, index_b) + rij(3)*fscalar
211 END IF
212
213 END DO bond
214
215 shell: DO ishell = 1, nshell
216 index_a = shell_list(ishell)%a + first_atom - 1
217 index_b = particle_set(index_a)%shell_index
218 rij = core_particle_set(index_b)%r - shell_particle_set(index_b)%r
219 rij = pbc(rij, cell)
220 k2_spring = shell_list(ishell)%shell_kind%k2_spring
221 k4_spring = shell_list(ishell)%shell_kind%k4_spring
222 r2 = dot_product(rij, rij)
223 energy = 0.5_dp*(k2_spring + k4_spring*r2/12.0_dp)*r2
224 fscalar = k2_spring + k4_spring*r2/6.0_dp
225 pot_shell = pot_shell + energy
226 IF (atener) THEN
227 ener_a(index_a) = ener_a(index_a) + energy
228 END IF
229 core_particle_set(index_b)%f(1) = core_particle_set(index_b)%f(1) - rij(1)*fscalar
230 core_particle_set(index_b)%f(2) = core_particle_set(index_b)%f(2) - rij(2)*fscalar
231 core_particle_set(index_b)%f(3) = core_particle_set(index_b)%f(3) - rij(3)*fscalar
232 shell_particle_set(index_b)%f(1) = shell_particle_set(index_b)%f(1) + rij(1)*fscalar
233 shell_particle_set(index_b)%f(2) = shell_particle_set(index_b)%f(2) + rij(2)*fscalar
234 shell_particle_set(index_b)%f(3) = shell_particle_set(index_b)%f(3) + rij(3)*fscalar
235 ! Compute the pressure tensor, if requested
236 IF (use_virial) THEN
237 k1 = -rij*fscalar
238 CALL get_pv_bond(k1, rij, pv_bond)
239 END IF
240 IF (atstress) THEN
241 CALL get_pv_bond(k1, rij, pv_a(:, :, index_a))
242 END IF
243 END DO shell
244
245 urey_bradley: DO ibend = 1, nub
246 index_a = ub_list(ibend)%a + first_atom - 1
247 index_b = ub_list(ibend)%c + first_atom - 1
248 rij = particle_set(index_a)%r - particle_set(index_b)%r
249 rij = pbc(rij, cell)
250 CALL force_bonds(ub_list(ibend)%ub_kind%id_type, rij, &
251 ub_list(ibend)%ub_kind%r0, &
252 ub_list(ibend)%ub_kind%k, 0.0_dp, energy, fscalar)
253 pot_urey_bradley = pot_urey_bradley + energy
254 IF (atener) THEN
255 ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
256 ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
257 END IF
258 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
259 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
260 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
261 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
262 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
263 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
264
265 ! computing the pressure tensor
266 k2 = -rij*fscalar
267 IF (use_virial) CALL get_pv_bond(k2, rij, pv_urey_bradley)
268 IF (atstress) THEN
269 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_a))
270 CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_b))
271 END IF
272
273 ! the contribution from the ub. ONLY FOR DEBUG
274 IF (PRESENT(f_ub)) THEN
275 f_ub(:, index_a) = f_ub(:, index_a) - rij*fscalar
276 f_ub(:, index_b) = f_ub(:, index_b) + rij*fscalar
277 END IF
278
279 END DO urey_bradley
280
281 bend: DO ibend = 1, nbends
282 index_a = bend_list(ibend)%a + first_atom - 1
283 index_b = bend_list(ibend)%b + first_atom - 1
284 index_c = bend_list(ibend)%c + first_atom - 1
285 b12 = particle_set(index_a)%r - particle_set(index_b)%r
286 b32 = particle_set(index_c)%r - particle_set(index_b)%r
287 b12 = pbc(b12, cell)
288 b32 = pbc(b32, cell)
289 d12 = sqrt(dot_product(b12, b12))
290 id12 = 1.0_dp/d12
291 d32 = sqrt(dot_product(b32, b32))
292 id32 = 1.0_dp/d32
293 dist = dot_product(b12, b32)
294 theta = (dist*id12*id32)
295 IF (theta < -1.0_dp) theta = -1.0_dp
296 IF (theta > +1.0_dp) theta = +1.0_dp
297 theta = acos(theta)
298 CALL force_bends(bend_list(ibend)%bend_kind%id_type, &
299 b12, b32, d12, d32, id12, id32, dist, theta, &
300 bend_list(ibend)%bend_kind%theta0, &
301 bend_list(ibend)%bend_kind%k, &
302 bend_list(ibend)%bend_kind%cb, &
303 bend_list(ibend)%bend_kind%r012, &
304 bend_list(ibend)%bend_kind%r032, &
305 bend_list(ibend)%bend_kind%kbs12, &
306 bend_list(ibend)%bend_kind%kbs32, &
307 bend_list(ibend)%bend_kind%kss, &
308 bend_list(ibend)%bend_kind%legendre, &
309 g1, g2, g3, energy, fscalar)
310 pot_bend = pot_bend + energy
311 IF (atener) THEN
312 ener_a(index_a) = ener_a(index_a) + energy/3._dp
313 ener_a(index_b) = ener_a(index_b) + energy/3._dp
314 ener_a(index_c) = ener_a(index_c) + energy/3._dp
315 END IF
316 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + g1(1)*fscalar
317 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + g1(2)*fscalar
318 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + g1(3)*fscalar
319 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + g2(1)*fscalar
320 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + g2(2)*fscalar
321 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + g2(3)*fscalar
322 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + g3(1)*fscalar
323 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + g3(2)*fscalar
324 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + g3(3)*fscalar
325
326 ! computing the pressure tensor
327 k1 = fscalar*g1
328 k3 = fscalar*g3
329 IF (use_virial) CALL get_pv_bend(k1, k3, b12, b32, pv_bend)
330 IF (atstress) THEN
331 k1 = fscalar*g1/3._dp
332 k3 = fscalar*g3/3._dp
333 CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_a))
334 CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_b))
335 CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_c))
336 END IF
337
338 ! the contribution from the bends. ONLY FOR DEBUG
339 IF (PRESENT(f_bend)) THEN
340 f_bend(:, index_a) = f_bend(:, index_a) + fscalar*g1
341 f_bend(:, index_b) = f_bend(:, index_b) + fscalar*g2
342 f_bend(:, index_c) = f_bend(:, index_c) + fscalar*g3
343 END IF
344 END DO bend
345
346 torsion: DO itorsion = 1, ntorsions
347 index_a = torsion_list(itorsion)%a + first_atom - 1
348 index_b = torsion_list(itorsion)%b + first_atom - 1
349 index_c = torsion_list(itorsion)%c + first_atom - 1
350 index_d = torsion_list(itorsion)%d + first_atom - 1
351 t12 = particle_set(index_a)%r - particle_set(index_b)%r
352 t32 = particle_set(index_c)%r - particle_set(index_b)%r
353 t34 = particle_set(index_c)%r - particle_set(index_d)%r
354 t43 = particle_set(index_d)%r - particle_set(index_c)%r
355 t12 = pbc(t12, cell)
356 t32 = pbc(t32, cell)
357 t34 = pbc(t34, cell)
358 t43 = pbc(t43, cell)
359 ! t12 x t32
360 tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
361 tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
362 tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
363 ! t32 x t34
364 tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
365 tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
366 tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
367 sm = sqrt(dot_product(tm, tm))
368 ism = 1.0_dp/sm
369 sn = sqrt(dot_product(tn, tn))
370 isn = 1.0_dp/sn
371 s32 = sqrt(dot_product(t32, t32))
372 is32 = 1.0_dp/s32
373 dist1 = dot_product(t12, t32)
374 dist2 = dot_product(t34, t32)
375 DO imul = 1, torsion_list(itorsion)%torsion_kind%nmul
376 CALL force_torsions(torsion_list(itorsion)%torsion_kind%id_type, &
377 s32, is32, ism, isn, dist1, dist2, tm, tn, t12, &
378 torsion_list(itorsion)%torsion_kind%k(imul), &
379 torsion_list(itorsion)%torsion_kind%phi0(imul), &
380 torsion_list(itorsion)%torsion_kind%m(imul), &
381 gt1, gt2, gt3, gt4, energy, fscalar)
382 pot_torsion = pot_torsion + energy
383 IF (atener) THEN
384 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
385 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
386 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
387 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
388 END IF
389 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
390 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
391 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
392 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
393 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
394 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
395 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
396 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
397 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
398 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
399 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
400 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
401
402 ! computing the pressure tensor
403 k1 = fscalar*gt1
404 k3 = fscalar*gt3
405 k4 = fscalar*gt4
406 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_torsion)
407 IF (atstress) THEN
408 k1 = 0.25_dp*fscalar*gt1
409 k3 = 0.25_dp*fscalar*gt3
410 k4 = 0.25_dp*fscalar*gt4
411 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a))
412 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b))
413 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c))
414 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_d))
415 END IF
416
417 ! the contribution from the torsions. ONLY FOR DEBUG
418 IF (PRESENT(f_torsion)) THEN
419 f_torsion(:, index_a) = f_torsion(:, index_a) + fscalar*gt1
420 f_torsion(:, index_b) = f_torsion(:, index_b) + fscalar*gt2
421 f_torsion(:, index_c) = f_torsion(:, index_c) + fscalar*gt3
422 f_torsion(:, index_d) = f_torsion(:, index_d) + fscalar*gt4
423 END IF
424 END DO
425 END DO torsion
426
427 imp_torsion: DO itorsion = 1, nimptors
428 index_a = impr_list(itorsion)%a + first_atom - 1
429 index_b = impr_list(itorsion)%b + first_atom - 1
430 index_c = impr_list(itorsion)%c + first_atom - 1
431 index_d = impr_list(itorsion)%d + first_atom - 1
432 t12 = particle_set(index_a)%r - particle_set(index_b)%r
433 t32 = particle_set(index_c)%r - particle_set(index_b)%r
434 t34 = particle_set(index_c)%r - particle_set(index_d)%r
435 t43 = particle_set(index_d)%r - particle_set(index_c)%r
436 t12 = pbc(t12, cell)
437 t32 = pbc(t32, cell)
438 t34 = pbc(t34, cell)
439 t43 = pbc(t43, cell)
440 ! t12 x t32
441 tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
442 tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
443 tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
444 ! t32 x t34
445 tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
446 tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
447 tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
448 sm = sqrt(dot_product(tm, tm))
449 ism = 1.0_dp/sm
450 sn = sqrt(dot_product(tn, tn))
451 isn = 1.0_dp/sn
452 s32 = sqrt(dot_product(t32, t32))
453 is32 = 1.0_dp/s32
454 dist1 = dot_product(t12, t32)
455 dist2 = dot_product(t34, t32)
456 CALL force_imp_torsions(impr_list(itorsion)%impr_kind%id_type, &
457 s32, is32, ism, isn, dist1, dist2, tm, tn, t12, &
458 impr_list(itorsion)%impr_kind%k, &
459 impr_list(itorsion)%impr_kind%phi0, &
460 gt1, gt2, gt3, gt4, energy, fscalar)
461 pot_imp_torsion = pot_imp_torsion + energy
462 IF (atener) THEN
463 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
464 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
465 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
466 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
467 END IF
468 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
469 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
470 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
471 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
472 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
473 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
474 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
475 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
476 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
477 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
478 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
479 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
480
481 ! computing the pressure tensor
482 k1 = fscalar*gt1
483 k3 = fscalar*gt3
484 k4 = fscalar*gt4
485 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_imp_torsion)
486 IF (atstress) THEN
487 k1 = 0.25_dp*fscalar*gt1
488 k3 = 0.25_dp*fscalar*gt3
489 k4 = 0.25_dp*fscalar*gt4
490 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a))
491 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b))
492 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c))
493 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_d))
494 END IF
495
496 ! the contribution from the torsions. ONLY FOR DEBUG
497 IF (PRESENT(f_imptor)) THEN
498 f_imptor(:, index_a) = f_imptor(:, index_a) + fscalar*gt1
499 f_imptor(:, index_b) = f_imptor(:, index_b) + fscalar*gt2
500 f_imptor(:, index_c) = f_imptor(:, index_c) + fscalar*gt3
501 f_imptor(:, index_d) = f_imptor(:, index_d) + fscalar*gt4
502 END IF
503 END DO imp_torsion
504
505 opbend: DO iopbend = 1, nopbends
506 index_a = opbend_list(iopbend)%a + first_atom - 1
507 index_b = opbend_list(iopbend)%b + first_atom - 1
508 index_c = opbend_list(iopbend)%c + first_atom - 1
509 index_d = opbend_list(iopbend)%d + first_atom - 1
510
511 t12 = particle_set(index_a)%r - particle_set(index_b)%r
512 t32 = particle_set(index_c)%r - particle_set(index_b)%r
513 t34 = particle_set(index_c)%r - particle_set(index_d)%r
514 t43 = particle_set(index_d)%r - particle_set(index_c)%r
515 t41 = particle_set(index_d)%r - particle_set(index_a)%r
516 t42 = pbc(t41 + t12, cell)
517 t12 = pbc(t12, cell)
518 t32 = pbc(t32, cell)
519 t41 = pbc(t41, cell)
520 t43 = pbc(t43, cell)
521 ! tm = t32 x t12
522 tm(1) = t32(2)*t12(3) - t12(2)*t32(3)
523 tm(2) = -t32(1)*t12(3) + t12(1)*t32(3)
524 tm(3) = t32(1)*t12(2) - t12(1)*t32(2)
525 sm = sqrt(dot_product(tm, tm))
526 s32 = sqrt(dot_product(t32, t32))
527 CALL force_opbends(opbend_list(iopbend)%opbend_kind%id_type, &
528 s32, tm, t41, t42, t43, &
529 opbend_list(iopbend)%opbend_kind%k, &
530 opbend_list(iopbend)%opbend_kind%phi0, &
531 gt1, gt2, gt3, gt4, energy, fscalar)
532 pot_opbend = pot_opbend + energy
533 IF (atener) THEN
534 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
535 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
536 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
537 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
538 END IF
539 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
540 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
541 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
542 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
543 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
544 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
545 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
546 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
547 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
548 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
549 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
550 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
551
552 ! computing the pressure tensor
553 k1 = fscalar*gt1
554 k3 = fscalar*gt3
555 k4 = fscalar*gt4
556
557 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_opbend)
558 IF (atstress) THEN
559 k1 = 0.25_dp*fscalar*gt1
560 k3 = 0.25_dp*fscalar*gt3
561 k4 = 0.25_dp*fscalar*gt4
562 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a))
563 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b))
564 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c))
565 CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_d))
566 END IF
567
568 ! the contribution from the opbends. ONLY FOR DEBUG
569 IF (PRESENT(f_opbend)) THEN
570 f_opbend(:, index_a) = f_opbend(:, index_a) + fscalar*gt1
571 f_opbend(:, index_b) = f_opbend(:, index_b) + fscalar*gt2
572 f_opbend(:, index_c) = f_opbend(:, index_c) + fscalar*gt3
573 f_opbend(:, index_d) = f_opbend(:, index_d) + fscalar*gt4
574 END IF
575 END DO opbend
576 END DO
577 END DO mol
578
579 CALL timestop(handle)
580
581 END SUBROUTINE force_intra_control
582
583END MODULE fist_intra_force
584
cell_types::pbc
Definition cell_types.F:92
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
fist_intra_force
Definition fist_intra_force.F:14
fist_intra_force::force_intra_control
subroutine, public force_intra_control(molecule_set, molecule_kind_set, local_molecules, particle_set, shell_particle_set, core_particle_set, pot_bond, pot_bend, pot_urey_bradley, pot_torsion, pot_imp_torsion, pot_opbend, pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, pv_imp_torsion, pv_opbend, f_bond, f_bend, f_torsion, f_ub, f_imptor, f_opbend, cell, use_virial, atprop_env)
Computes the the intramolecular energies, forces, and pressure tensors.
Definition fist_intra_force.F:87
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
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:277
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
molecule_kind_types::bend_type
Definition molecule_kind_types.F:76
molecule_kind_types::bond_type
Definition molecule_kind_types.F:69
molecule_kind_types::impr_type
Definition molecule_kind_types.F:97
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_kind_types::opbend_type
Definition molecule_kind_types.F:104
molecule_kind_types::shell_type
Definition molecule_kind_types.F:62
molecule_kind_types::torsion_type
Definition molecule_kind_types.F:90
molecule_kind_types::ub_type
Definition molecule_kind_types.F:83
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35