(git:ed6f26b)
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-2025 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
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 TYPE(bend_type), POINTER :: bend_list(:)
120 TYPE(bond_type), POINTER :: bond_list(:)
121 TYPE(cp_logger_type), POINTER :: logger
122 TYPE(impr_type), POINTER :: impr_list(:)
123 TYPE(molecule_kind_type), POINTER :: molecule_kind
124 TYPE(molecule_type), POINTER :: molecule
125 TYPE(opbend_type), POINTER :: opbend_list(:)
126 TYPE(shell_type), POINTER :: shell_list(:)
127 TYPE(torsion_type), POINTER :: torsion_list(:)
128 TYPE(ub_type), POINTER :: ub_list(:)
129
130 CALL timeset(routinen, handle)
131 NULLIFY (logger)
132 logger => cp_get_default_logger()
133
134 IF (PRESENT(f_bond)) f_bond = 0.0_dp
135 IF (PRESENT(f_bend)) f_bend = 0.0_dp
136 IF (PRESENT(f_torsion)) f_torsion = 0.0_dp
137 IF (PRESENT(f_imptor)) f_imptor = 0.0_dp
138 IF (PRESENT(f_ub)) f_ub = 0.0_dp
139
140 pot_bond = 0.0_dp
141 pot_bend = 0.0_dp
142 pot_urey_bradley = 0.0_dp
143 pot_torsion = 0.0_dp
144 pot_imp_torsion = 0.0_dp
145 pot_opbend = 0.0_dp
146 pot_shell = 0.0_dp
147
148 atener = atprop_env%energy
149 IF (atener) ener_a => atprop_env%atener
150
151 nkind = SIZE(molecule_kind_set)
152 mol: DO ikind = 1, nkind
153 nmol_per_kind = local_molecules%n_el(ikind)
154
155 DO imol = 1, nmol_per_kind
156 i = local_molecules%list(ikind)%array(imol)
157 molecule => molecule_set(i)
158 molecule_kind => molecule%molecule_kind
159 CALL get_molecule_kind(molecule_kind, nbend=nbends, nbond=nbonds, &
160 nimpr=nimptors, nub=nub, ntorsion=ntorsions, &
161 nopbend=nopbends, nshell=nshell, &
162 bond_list=bond_list, ub_list=ub_list, &
163 bend_list=bend_list, torsion_list=torsion_list, &
164 impr_list=impr_list, opbend_list=opbend_list, &
165 shell_list=shell_list)
166
167 CALL get_molecule(molecule, first_atom=first_atom)
168
169 bond: DO ibond = 1, nbonds
170 index_a = bond_list(ibond)%a + first_atom - 1
171 index_b = bond_list(ibond)%b + first_atom - 1
172 rij = particle_set(index_a)%r - particle_set(index_b)%r
173 rij = pbc(rij, cell)
174 CALL force_bonds(bond_list(ibond)%bond_kind%id_type, rij, &
175 bond_list(ibond)%bond_kind%r0, &
176 bond_list(ibond)%bond_kind%k, &
177 bond_list(ibond)%bond_kind%cs, &
178 energy, fscalar)
179 pot_bond = pot_bond + energy
180 IF (atener) THEN
181 ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
182 ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
183 END IF
184
185 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
186 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
187 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
188 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
189 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
190 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
191
192 ! computing the pressure tensor
193 k2 = -rij*fscalar
194 IF (use_virial) CALL get_pv_bond(k2, rij, pv_bond)
195
196 ! the contribution from the bonds. ONLY FOR DEBUG
197 IF (PRESENT(f_bond)) THEN
198 f_bond(1, index_a) = f_bond(1, index_a) - rij(1)*fscalar
199 f_bond(2, index_a) = f_bond(2, index_a) - rij(2)*fscalar
200 f_bond(3, index_a) = f_bond(3, index_a) - rij(3)*fscalar
201 f_bond(1, index_b) = f_bond(1, index_b) + rij(1)*fscalar
202 f_bond(2, index_b) = f_bond(2, index_b) + rij(2)*fscalar
203 f_bond(3, index_b) = f_bond(3, index_b) + rij(3)*fscalar
204 END IF
205
206 END DO bond
207
208 shell: DO ishell = 1, nshell
209 index_a = shell_list(ishell)%a + first_atom - 1
210 index_b = particle_set(index_a)%shell_index
211 rij = core_particle_set(index_b)%r - shell_particle_set(index_b)%r
212 rij = pbc(rij, cell)
213 k2_spring = shell_list(ishell)%shell_kind%k2_spring
214 k4_spring = shell_list(ishell)%shell_kind%k4_spring
215 r2 = dot_product(rij, rij)
216 energy = 0.5_dp*(k2_spring + k4_spring*r2/12.0_dp)*r2
217 fscalar = k2_spring + k4_spring*r2/6.0_dp
218 pot_shell = pot_shell + energy
219 IF (atener) THEN
220 ener_a(index_a) = ener_a(index_a) + energy
221 END IF
222 core_particle_set(index_b)%f(1) = core_particle_set(index_b)%f(1) - rij(1)*fscalar
223 core_particle_set(index_b)%f(2) = core_particle_set(index_b)%f(2) - rij(2)*fscalar
224 core_particle_set(index_b)%f(3) = core_particle_set(index_b)%f(3) - rij(3)*fscalar
225 shell_particle_set(index_b)%f(1) = shell_particle_set(index_b)%f(1) + rij(1)*fscalar
226 shell_particle_set(index_b)%f(2) = shell_particle_set(index_b)%f(2) + rij(2)*fscalar
227 shell_particle_set(index_b)%f(3) = shell_particle_set(index_b)%f(3) + rij(3)*fscalar
228 ! Compute the pressure tensor, if requested
229 IF (use_virial) THEN
230 k1 = -rij*fscalar
231 CALL get_pv_bond(k1, rij, pv_bond)
232 END IF
233 END DO shell
234
235 urey_bradley: DO ibend = 1, nub
236 index_a = ub_list(ibend)%a + first_atom - 1
237 index_b = ub_list(ibend)%c + first_atom - 1
238 rij = particle_set(index_a)%r - particle_set(index_b)%r
239 rij = pbc(rij, cell)
240 CALL force_bonds(ub_list(ibend)%ub_kind%id_type, rij, &
241 ub_list(ibend)%ub_kind%r0, &
242 ub_list(ibend)%ub_kind%k, 0.0_dp, energy, fscalar)
243 pot_urey_bradley = pot_urey_bradley + energy
244 IF (atener) THEN
245 ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
246 ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
247 END IF
248 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
249 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
250 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
251 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
252 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
253 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
254
255 ! computing the pressure tensor
256 k2 = -rij*fscalar
257 IF (use_virial) CALL get_pv_bond(k2, rij, pv_urey_bradley)
258
259 ! the contribution from the ub. ONLY FOR DEBUG
260 IF (PRESENT(f_ub)) THEN
261 f_ub(:, index_a) = f_ub(:, index_a) - rij*fscalar
262 f_ub(:, index_b) = f_ub(:, index_b) + rij*fscalar
263 END IF
264
265 END DO urey_bradley
266
267 bend: DO ibend = 1, nbends
268 index_a = bend_list(ibend)%a + first_atom - 1
269 index_b = bend_list(ibend)%b + first_atom - 1
270 index_c = bend_list(ibend)%c + first_atom - 1
271 b12 = particle_set(index_a)%r - particle_set(index_b)%r
272 b32 = particle_set(index_c)%r - particle_set(index_b)%r
273 b12 = pbc(b12, cell)
274 b32 = pbc(b32, cell)
275 d12 = sqrt(dot_product(b12, b12))
276 id12 = 1.0_dp/d12
277 d32 = sqrt(dot_product(b32, b32))
278 id32 = 1.0_dp/d32
279 dist = dot_product(b12, b32)
280 theta = (dist*id12*id32)
281 IF (theta < -1.0_dp) theta = -1.0_dp
282 IF (theta > +1.0_dp) theta = +1.0_dp
283 theta = acos(theta)
284 CALL force_bends(bend_list(ibend)%bend_kind%id_type, &
285 b12, b32, d12, d32, id12, id32, dist, theta, &
286 bend_list(ibend)%bend_kind%theta0, &
287 bend_list(ibend)%bend_kind%k, &
288 bend_list(ibend)%bend_kind%cb, &
289 bend_list(ibend)%bend_kind%r012, &
290 bend_list(ibend)%bend_kind%r032, &
291 bend_list(ibend)%bend_kind%kbs12, &
292 bend_list(ibend)%bend_kind%kbs32, &
293 bend_list(ibend)%bend_kind%kss, &
294 bend_list(ibend)%bend_kind%legendre, &
295 g1, g2, g3, energy, fscalar)
296 pot_bend = pot_bend + energy
297 IF (atener) THEN
298 ener_a(index_a) = ener_a(index_a) + energy/3._dp
299 ener_a(index_b) = ener_a(index_b) + energy/3._dp
300 ener_a(index_c) = ener_a(index_c) + energy/3._dp
301 END IF
302 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + g1(1)*fscalar
303 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + g1(2)*fscalar
304 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + g1(3)*fscalar
305 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + g2(1)*fscalar
306 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + g2(2)*fscalar
307 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + g2(3)*fscalar
308 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + g3(1)*fscalar
309 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + g3(2)*fscalar
310 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + g3(3)*fscalar
311
312 ! computing the pressure tensor
313 k1 = fscalar*g1
314 k3 = fscalar*g3
315 IF (use_virial) CALL get_pv_bend(k1, k3, b12, b32, pv_bend)
316
317 ! the contribution from the bends. ONLY FOR DEBUG
318 IF (PRESENT(f_bend)) THEN
319 f_bend(:, index_a) = f_bend(:, index_a) + fscalar*g1
320 f_bend(:, index_b) = f_bend(:, index_b) + fscalar*g2
321 f_bend(:, index_c) = f_bend(:, index_c) + fscalar*g3
322 END IF
323 END DO bend
324
325 torsion: DO itorsion = 1, ntorsions
326 index_a = torsion_list(itorsion)%a + first_atom - 1
327 index_b = torsion_list(itorsion)%b + first_atom - 1
328 index_c = torsion_list(itorsion)%c + first_atom - 1
329 index_d = torsion_list(itorsion)%d + first_atom - 1
330 t12 = particle_set(index_a)%r - particle_set(index_b)%r
331 t32 = particle_set(index_c)%r - particle_set(index_b)%r
332 t34 = particle_set(index_c)%r - particle_set(index_d)%r
333 t43 = particle_set(index_d)%r - particle_set(index_c)%r
334 t12 = pbc(t12, cell)
335 t32 = pbc(t32, cell)
336 t34 = pbc(t34, cell)
337 t43 = pbc(t43, cell)
338 ! t12 x t32
339 tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
340 tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
341 tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
342 ! t32 x t34
343 tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
344 tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
345 tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
346 sm = sqrt(dot_product(tm, tm))
347 ism = 1.0_dp/sm
348 sn = sqrt(dot_product(tn, tn))
349 isn = 1.0_dp/sn
350 s32 = sqrt(dot_product(t32, t32))
351 is32 = 1.0_dp/s32
352 dist1 = dot_product(t12, t32)
353 dist2 = dot_product(t34, t32)
354 DO imul = 1, torsion_list(itorsion)%torsion_kind%nmul
355 CALL force_torsions(torsion_list(itorsion)%torsion_kind%id_type, &
356 s32, is32, ism, isn, dist1, dist2, tm, tn, t12, &
357 torsion_list(itorsion)%torsion_kind%k(imul), &
358 torsion_list(itorsion)%torsion_kind%phi0(imul), &
359 torsion_list(itorsion)%torsion_kind%m(imul), &
360 gt1, gt2, gt3, gt4, energy, fscalar)
361 pot_torsion = pot_torsion + energy
362 IF (atener) THEN
363 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
364 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
365 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
366 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
367 END IF
368 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
369 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
370 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
371 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
372 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
373 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
374 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
375 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
376 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
377 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
378 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
379 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
380
381 ! computing the pressure tensor
382 k1 = fscalar*gt1
383 k3 = fscalar*gt3
384 k4 = fscalar*gt4
385 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_torsion)
386
387 ! the contribution from the torsions. ONLY FOR DEBUG
388 IF (PRESENT(f_torsion)) THEN
389 f_torsion(:, index_a) = f_torsion(:, index_a) + fscalar*gt1
390 f_torsion(:, index_b) = f_torsion(:, index_b) + fscalar*gt2
391 f_torsion(:, index_c) = f_torsion(:, index_c) + fscalar*gt3
392 f_torsion(:, index_d) = f_torsion(:, index_d) + fscalar*gt4
393 END IF
394 END DO
395 END DO torsion
396
397 imp_torsion: DO itorsion = 1, nimptors
398 index_a = impr_list(itorsion)%a + first_atom - 1
399 index_b = impr_list(itorsion)%b + first_atom - 1
400 index_c = impr_list(itorsion)%c + first_atom - 1
401 index_d = impr_list(itorsion)%d + first_atom - 1
402 t12 = particle_set(index_a)%r - particle_set(index_b)%r
403 t32 = particle_set(index_c)%r - particle_set(index_b)%r
404 t34 = particle_set(index_c)%r - particle_set(index_d)%r
405 t43 = particle_set(index_d)%r - particle_set(index_c)%r
406 t12 = pbc(t12, cell)
407 t32 = pbc(t32, cell)
408 t34 = pbc(t34, cell)
409 t43 = pbc(t43, cell)
410 ! t12 x t32
411 tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
412 tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
413 tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
414 ! t32 x t34
415 tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
416 tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
417 tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
418 sm = sqrt(dot_product(tm, tm))
419 ism = 1.0_dp/sm
420 sn = sqrt(dot_product(tn, tn))
421 isn = 1.0_dp/sn
422 s32 = sqrt(dot_product(t32, t32))
423 is32 = 1.0_dp/s32
424 dist1 = dot_product(t12, t32)
425 dist2 = dot_product(t34, t32)
426 CALL force_imp_torsions(impr_list(itorsion)%impr_kind%id_type, &
427 s32, is32, ism, isn, dist1, dist2, tm, tn, t12, &
428 impr_list(itorsion)%impr_kind%k, &
429 impr_list(itorsion)%impr_kind%phi0, &
430 gt1, gt2, gt3, gt4, energy, fscalar)
431 pot_imp_torsion = pot_imp_torsion + energy
432 IF (atener) THEN
433 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
434 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
435 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
436 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
437 END IF
438 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
439 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
440 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
441 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
442 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
443 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
444 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
445 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
446 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
447 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
448 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
449 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
450
451 ! computing the pressure tensor
452 k1 = fscalar*gt1
453 k3 = fscalar*gt3
454 k4 = fscalar*gt4
455 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_imp_torsion)
456
457 ! the contribution from the torsions. ONLY FOR DEBUG
458 IF (PRESENT(f_imptor)) THEN
459 f_imptor(:, index_a) = f_imptor(:, index_a) + fscalar*gt1
460 f_imptor(:, index_b) = f_imptor(:, index_b) + fscalar*gt2
461 f_imptor(:, index_c) = f_imptor(:, index_c) + fscalar*gt3
462 f_imptor(:, index_d) = f_imptor(:, index_d) + fscalar*gt4
463 END IF
464 END DO imp_torsion
465
466 opbend: DO iopbend = 1, nopbends
467 index_a = opbend_list(iopbend)%a + first_atom - 1
468 index_b = opbend_list(iopbend)%b + first_atom - 1
469 index_c = opbend_list(iopbend)%c + first_atom - 1
470 index_d = opbend_list(iopbend)%d + first_atom - 1
471
472 t12 = particle_set(index_a)%r - particle_set(index_b)%r
473 t32 = particle_set(index_c)%r - particle_set(index_b)%r
474 t34 = particle_set(index_c)%r - particle_set(index_d)%r
475 t43 = particle_set(index_d)%r - particle_set(index_c)%r
476 t41 = particle_set(index_d)%r - particle_set(index_a)%r
477 t42 = pbc(t41 + t12, cell)
478 t12 = pbc(t12, cell)
479 t32 = pbc(t32, cell)
480 t41 = pbc(t41, cell)
481 t43 = pbc(t43, cell)
482 ! tm = t32 x t12
483 tm(1) = t32(2)*t12(3) - t12(2)*t32(3)
484 tm(2) = -t32(1)*t12(3) + t12(1)*t32(3)
485 tm(3) = t32(1)*t12(2) - t12(1)*t32(2)
486 sm = sqrt(dot_product(tm, tm))
487 s32 = sqrt(dot_product(t32, t32))
488 CALL force_opbends(opbend_list(iopbend)%opbend_kind%id_type, &
489 s32, tm, t41, t42, t43, &
490 opbend_list(iopbend)%opbend_kind%k, &
491 opbend_list(iopbend)%opbend_kind%phi0, &
492 gt1, gt2, gt3, gt4, energy, fscalar)
493 pot_opbend = pot_opbend + energy
494 IF (atener) THEN
495 ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
496 ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
497 ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
498 ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
499 END IF
500 particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
501 particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
502 particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
503 particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
504 particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
505 particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
506 particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
507 particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
508 particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
509 particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
510 particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
511 particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
512
513 ! computing the pressure tensor
514 k1 = fscalar*gt1
515 k3 = fscalar*gt3
516 k4 = fscalar*gt4
517
518 IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_opbend)
519
520 ! the contribution from the opbends. ONLY FOR DEBUG
521 IF (PRESENT(f_opbend)) THEN
522 f_opbend(:, index_a) = f_opbend(:, index_a) + fscalar*gt1
523 f_opbend(:, index_b) = f_opbend(:, index_b) + fscalar*gt2
524 f_opbend(:, index_c) = f_opbend(:, index_c) + fscalar*gt3
525 f_opbend(:, index_d) = f_opbend(:, index_d) + fscalar*gt4
526 END IF
527 END DO opbend
528 END DO
529 END DO mol
530
531 CALL timestop(handle)
532
533 END SUBROUTINE force_intra_control
534
535END MODULE fist_intra_force
536
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