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 : !> Variable names changed (DG) 05-Dec-2000
12 : !> \author CJM
13 : ! **************************************************************************************************
14 : MODULE 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 :
45 : CONTAINS
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 166176 : 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 83088 : pot_opbend, pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, &
85 249264 : pv_imp_torsion, pv_opbend, f_bond, f_bend, f_torsion, f_ub, &
86 166176 : 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 83088 : REAL(KIND=dp), DIMENSION(:), POINTER :: ener_a
119 83088 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pv_a
120 83088 : TYPE(bend_type), POINTER :: bend_list(:)
121 83088 : TYPE(bond_type), POINTER :: bond_list(:)
122 : TYPE(cp_logger_type), POINTER :: logger
123 83088 : TYPE(impr_type), POINTER :: impr_list(:)
124 : TYPE(molecule_kind_type), POINTER :: molecule_kind
125 : TYPE(molecule_type), POINTER :: molecule
126 83088 : TYPE(opbend_type), POINTER :: opbend_list(:)
127 83088 : TYPE(shell_type), POINTER :: shell_list(:)
128 83088 : TYPE(torsion_type), POINTER :: torsion_list(:)
129 83088 : TYPE(ub_type), POINTER :: ub_list(:)
130 :
131 83088 : CALL timeset(routineN, handle)
132 83088 : NULLIFY (logger)
133 83088 : logger => cp_get_default_logger()
134 :
135 83088 : IF (PRESENT(f_bond)) f_bond = 0.0_dp
136 83088 : IF (PRESENT(f_bend)) f_bend = 0.0_dp
137 83088 : IF (PRESENT(f_torsion)) f_torsion = 0.0_dp
138 83088 : IF (PRESENT(f_imptor)) f_imptor = 0.0_dp
139 83088 : IF (PRESENT(f_ub)) f_ub = 0.0_dp
140 :
141 83088 : pot_bond = 0.0_dp
142 83088 : pot_bend = 0.0_dp
143 83088 : pot_urey_bradley = 0.0_dp
144 83088 : pot_torsion = 0.0_dp
145 83088 : pot_imp_torsion = 0.0_dp
146 83088 : pot_opbend = 0.0_dp
147 83088 : pot_shell = 0.0_dp
148 :
149 83088 : atener = atprop_env%energy
150 83088 : IF (atener) ener_a => atprop_env%atener
151 83088 : atstress = atprop_env%stress
152 83088 : IF (atstress) pv_a => atprop_env%atstress
153 :
154 83088 : nkind = SIZE(molecule_kind_set)
155 1407323 : MOL: DO ikind = 1, nkind
156 1324235 : nmol_per_kind = local_molecules%n_el(ikind)
157 :
158 3457916 : DO imol = 1, nmol_per_kind
159 2050593 : i = local_molecules%list(ikind)%array(imol)
160 2050593 : molecule => molecule_set(i)
161 2050593 : 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 2050593 : shell_list=shell_list)
169 :
170 2050593 : CALL get_molecule(molecule, first_atom=first_atom)
171 :
172 4820693 : BOND: DO ibond = 1, nbonds
173 2770100 : index_a = bond_list(ibond)%a + first_atom - 1
174 2770100 : index_b = bond_list(ibond)%b + first_atom - 1
175 11080400 : rij = particle_set(index_a)%r - particle_set(index_b)%r
176 11080400 : 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 2770100 : energy, fscalar)
182 2770100 : pot_bond = pot_bond + energy
183 2770100 : IF (atener) THEN
184 606 : ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
185 606 : ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
186 : END IF
187 :
188 2770100 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
189 2770100 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
190 2770100 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
191 2770100 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
192 2770100 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
193 2770100 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
194 :
195 : ! computing the pressure tensor
196 11080400 : k2 = -rij*fscalar
197 2770100 : IF (use_virial) CALL get_pv_bond(k2, rij, pv_bond)
198 2770100 : IF (atstress) THEN
199 2424 : CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_a))
200 2424 : 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 4820693 : IF (PRESENT(f_bond)) THEN
205 0 : f_bond(1, index_a) = f_bond(1, index_a) - rij(1)*fscalar
206 0 : f_bond(2, index_a) = f_bond(2, index_a) - rij(2)*fscalar
207 0 : f_bond(3, index_a) = f_bond(3, index_a) - rij(3)*fscalar
208 0 : f_bond(1, index_b) = f_bond(1, index_b) + rij(1)*fscalar
209 0 : f_bond(2, index_b) = f_bond(2, index_b) + rij(2)*fscalar
210 0 : f_bond(3, index_b) = f_bond(3, index_b) + rij(3)*fscalar
211 : END IF
212 :
213 : END DO BOND
214 :
215 2466480 : SHELL: DO ishell = 1, nshell
216 415887 : index_a = shell_list(ishell)%a + first_atom - 1
217 415887 : index_b = particle_set(index_a)%shell_index
218 1663548 : rij = core_particle_set(index_b)%r - shell_particle_set(index_b)%r
219 1663548 : rij = pbc(rij, cell)
220 415887 : k2_spring = shell_list(ishell)%shell_kind%k2_spring
221 415887 : k4_spring = shell_list(ishell)%shell_kind%k4_spring
222 1663548 : r2 = DOT_PRODUCT(rij, rij)
223 415887 : energy = 0.5_dp*(k2_spring + k4_spring*r2/12.0_dp)*r2
224 415887 : fscalar = k2_spring + k4_spring*r2/6.0_dp
225 415887 : pot_shell = pot_shell + energy
226 415887 : IF (atener) THEN
227 1080 : ener_a(index_a) = ener_a(index_a) + energy
228 : END IF
229 415887 : core_particle_set(index_b)%f(1) = core_particle_set(index_b)%f(1) - rij(1)*fscalar
230 415887 : core_particle_set(index_b)%f(2) = core_particle_set(index_b)%f(2) - rij(2)*fscalar
231 415887 : core_particle_set(index_b)%f(3) = core_particle_set(index_b)%f(3) - rij(3)*fscalar
232 415887 : shell_particle_set(index_b)%f(1) = shell_particle_set(index_b)%f(1) + rij(1)*fscalar
233 415887 : shell_particle_set(index_b)%f(2) = shell_particle_set(index_b)%f(2) + rij(2)*fscalar
234 415887 : 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 415887 : IF (use_virial) THEN
237 1366840 : k1 = -rij*fscalar
238 341710 : CALL get_pv_bond(k1, rij, pv_bond)
239 : END IF
240 2466480 : IF (atstress) THEN
241 0 : CALL get_pv_bond(k1, rij, pv_a(:, :, index_a))
242 : END IF
243 : END DO SHELL
244 :
245 2158796 : UREY_BRADLEY: DO ibend = 1, nub
246 108203 : index_a = ub_list(ibend)%a + first_atom - 1
247 108203 : index_b = ub_list(ibend)%c + first_atom - 1
248 432812 : rij = particle_set(index_a)%r - particle_set(index_b)%r
249 432812 : rij = pbc(rij, cell)
250 : CALL force_bonds(ub_list(ibend)%ub_kind%id_type, rij, &
251 : ub_list(ibend)%ub_kind%r0, &
252 108203 : ub_list(ibend)%ub_kind%k, 0.0_dp, energy, fscalar)
253 108203 : pot_urey_bradley = pot_urey_bradley + energy
254 108203 : IF (atener) THEN
255 0 : ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
256 0 : ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
257 : END IF
258 108203 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
259 108203 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
260 108203 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
261 108203 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
262 108203 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
263 108203 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
264 :
265 : ! computing the pressure tensor
266 432812 : k2 = -rij*fscalar
267 108203 : IF (use_virial) CALL get_pv_bond(k2, rij, pv_urey_bradley)
268 108203 : IF (atstress) THEN
269 0 : CALL get_pv_bond(0.5_dp*k2, rij, pv_a(:, :, index_a))
270 0 : 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 2158796 : IF (PRESENT(f_ub)) THEN
275 0 : f_ub(:, index_a) = f_ub(:, index_a) - rij*fscalar
276 0 : f_ub(:, index_b) = f_ub(:, index_b) + rij*fscalar
277 : END IF
278 :
279 : END DO UREY_BRADLEY
280 :
281 3400251 : BEND: DO ibend = 1, nbends
282 1349658 : index_a = bend_list(ibend)%a + first_atom - 1
283 1349658 : index_b = bend_list(ibend)%b + first_atom - 1
284 1349658 : index_c = bend_list(ibend)%c + first_atom - 1
285 5398632 : b12 = particle_set(index_a)%r - particle_set(index_b)%r
286 5398632 : b32 = particle_set(index_c)%r - particle_set(index_b)%r
287 5398632 : b12 = pbc(b12, cell)
288 5398632 : b32 = pbc(b32, cell)
289 5398632 : d12 = SQRT(DOT_PRODUCT(b12, b12))
290 1349658 : id12 = 1.0_dp/d12
291 5398632 : d32 = SQRT(DOT_PRODUCT(b32, b32))
292 1349658 : id32 = 1.0_dp/d32
293 5398632 : dist = DOT_PRODUCT(b12, b32)
294 1349658 : theta = (dist*id12*id32)
295 1349658 : IF (theta < -1.0_dp) theta = -1.0_dp
296 1349658 : IF (theta > +1.0_dp) theta = +1.0_dp
297 1349658 : 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 1349658 : g1, g2, g3, energy, fscalar)
310 1349658 : pot_bend = pot_bend + energy
311 1349658 : IF (atener) THEN
312 303 : ener_a(index_a) = ener_a(index_a) + energy/3._dp
313 303 : ener_a(index_b) = ener_a(index_b) + energy/3._dp
314 303 : ener_a(index_c) = ener_a(index_c) + energy/3._dp
315 : END IF
316 1349658 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + g1(1)*fscalar
317 1349658 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + g1(2)*fscalar
318 1349658 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + g1(3)*fscalar
319 1349658 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + g2(1)*fscalar
320 1349658 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + g2(2)*fscalar
321 1349658 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + g2(3)*fscalar
322 1349658 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + g3(1)*fscalar
323 1349658 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + g3(2)*fscalar
324 1349658 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + g3(3)*fscalar
325 :
326 : ! computing the pressure tensor
327 5398632 : k1 = fscalar*g1
328 5398632 : k3 = fscalar*g3
329 1349658 : IF (use_virial) CALL get_pv_bend(k1, k3, b12, b32, pv_bend)
330 1349658 : IF (atstress) THEN
331 1212 : k1 = fscalar*g1/3._dp
332 1212 : k3 = fscalar*g3/3._dp
333 303 : CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_a))
334 303 : CALL get_pv_bend(k1, k3, b12, b32, pv_a(:, :, index_b))
335 303 : 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 3400251 : IF (PRESENT(f_bend)) THEN
340 0 : f_bend(:, index_a) = f_bend(:, index_a) + fscalar*g1
341 0 : f_bend(:, index_b) = f_bend(:, index_b) + fscalar*g2
342 0 : f_bend(:, index_c) = f_bend(:, index_c) + fscalar*g3
343 : END IF
344 : END DO BEND
345 :
346 2737899 : TORSION: DO itorsion = 1, ntorsions
347 687306 : index_a = torsion_list(itorsion)%a + first_atom - 1
348 687306 : index_b = torsion_list(itorsion)%b + first_atom - 1
349 687306 : index_c = torsion_list(itorsion)%c + first_atom - 1
350 687306 : index_d = torsion_list(itorsion)%d + first_atom - 1
351 2749224 : t12 = particle_set(index_a)%r - particle_set(index_b)%r
352 2749224 : t32 = particle_set(index_c)%r - particle_set(index_b)%r
353 2749224 : t34 = particle_set(index_c)%r - particle_set(index_d)%r
354 2749224 : t43 = particle_set(index_d)%r - particle_set(index_c)%r
355 2749224 : t12 = pbc(t12, cell)
356 2749224 : t32 = pbc(t32, cell)
357 2749224 : t34 = pbc(t34, cell)
358 2749224 : t43 = pbc(t43, cell)
359 : ! t12 x t32
360 687306 : tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
361 687306 : tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
362 687306 : tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
363 : ! t32 x t34
364 687306 : tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
365 687306 : tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
366 687306 : tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
367 2749224 : sm = SQRT(DOT_PRODUCT(tm, tm))
368 687306 : ism = 1.0_dp/sm
369 2749224 : sn = SQRT(DOT_PRODUCT(tn, tn))
370 687306 : isn = 1.0_dp/sn
371 2749224 : s32 = SQRT(DOT_PRODUCT(t32, t32))
372 687306 : is32 = 1.0_dp/s32
373 2749224 : dist1 = DOT_PRODUCT(t12, t32)
374 2749224 : dist2 = DOT_PRODUCT(t34, t32)
375 3483922 : 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 746023 : gt1, gt2, gt3, gt4, energy, fscalar)
382 746023 : pot_torsion = pot_torsion + energy
383 746023 : IF (atener) THEN
384 0 : ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
385 0 : ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
386 0 : ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
387 0 : ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
388 : END IF
389 746023 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
390 746023 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
391 746023 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
392 746023 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
393 746023 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
394 746023 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
395 746023 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
396 746023 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
397 746023 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
398 746023 : particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
399 746023 : particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
400 746023 : particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
401 :
402 : ! computing the pressure tensor
403 2984092 : k1 = fscalar*gt1
404 2984092 : k3 = fscalar*gt3
405 2984092 : k4 = fscalar*gt4
406 746023 : IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_torsion)
407 746023 : IF (atstress) THEN
408 0 : k1 = 0.25_dp*fscalar*gt1
409 0 : k3 = 0.25_dp*fscalar*gt3
410 0 : k4 = 0.25_dp*fscalar*gt4
411 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a))
412 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b))
413 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c))
414 0 : 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 1433329 : IF (PRESENT(f_torsion)) THEN
419 0 : f_torsion(:, index_a) = f_torsion(:, index_a) + fscalar*gt1
420 0 : f_torsion(:, index_b) = f_torsion(:, index_b) + fscalar*gt2
421 0 : f_torsion(:, index_c) = f_torsion(:, index_c) + fscalar*gt3
422 0 : f_torsion(:, index_d) = f_torsion(:, index_d) + fscalar*gt4
423 : END IF
424 : END DO
425 : END DO TORSION
426 :
427 2069900 : IMP_TORSION: DO itorsion = 1, nimptors
428 19307 : index_a = impr_list(itorsion)%a + first_atom - 1
429 19307 : index_b = impr_list(itorsion)%b + first_atom - 1
430 19307 : index_c = impr_list(itorsion)%c + first_atom - 1
431 19307 : index_d = impr_list(itorsion)%d + first_atom - 1
432 77228 : t12 = particle_set(index_a)%r - particle_set(index_b)%r
433 77228 : t32 = particle_set(index_c)%r - particle_set(index_b)%r
434 77228 : t34 = particle_set(index_c)%r - particle_set(index_d)%r
435 77228 : t43 = particle_set(index_d)%r - particle_set(index_c)%r
436 77228 : t12 = pbc(t12, cell)
437 77228 : t32 = pbc(t32, cell)
438 77228 : t34 = pbc(t34, cell)
439 77228 : t43 = pbc(t43, cell)
440 : ! t12 x t32
441 19307 : tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
442 19307 : tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
443 19307 : tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
444 : ! t32 x t34
445 19307 : tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
446 19307 : tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
447 19307 : tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
448 77228 : sm = SQRT(DOT_PRODUCT(tm, tm))
449 19307 : ism = 1.0_dp/sm
450 77228 : sn = SQRT(DOT_PRODUCT(tn, tn))
451 19307 : isn = 1.0_dp/sn
452 77228 : s32 = SQRT(DOT_PRODUCT(t32, t32))
453 19307 : is32 = 1.0_dp/s32
454 77228 : dist1 = DOT_PRODUCT(t12, t32)
455 77228 : 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 19307 : gt1, gt2, gt3, gt4, energy, fscalar)
461 19307 : pot_imp_torsion = pot_imp_torsion + energy
462 19307 : IF (atener) THEN
463 0 : ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
464 0 : ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
465 0 : ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
466 0 : ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
467 : END IF
468 19307 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
469 19307 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
470 19307 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
471 19307 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
472 19307 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
473 19307 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
474 19307 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
475 19307 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
476 19307 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
477 19307 : particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
478 19307 : particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
479 19307 : particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
480 :
481 : ! computing the pressure tensor
482 77228 : k1 = fscalar*gt1
483 77228 : k3 = fscalar*gt3
484 77228 : k4 = fscalar*gt4
485 19307 : IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_imp_torsion)
486 19307 : IF (atstress) THEN
487 0 : k1 = 0.25_dp*fscalar*gt1
488 0 : k3 = 0.25_dp*fscalar*gt3
489 0 : k4 = 0.25_dp*fscalar*gt4
490 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a))
491 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b))
492 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c))
493 0 : 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 2069900 : IF (PRESENT(f_imptor)) THEN
498 0 : f_imptor(:, index_a) = f_imptor(:, index_a) + fscalar*gt1
499 0 : f_imptor(:, index_b) = f_imptor(:, index_b) + fscalar*gt2
500 0 : f_imptor(:, index_c) = f_imptor(:, index_c) + fscalar*gt3
501 0 : f_imptor(:, index_d) = f_imptor(:, index_d) + fscalar*gt4
502 : END IF
503 : END DO IMP_TORSION
504 :
505 5425446 : OPBEND: DO iopbend = 1, nopbends
506 25 : index_a = opbend_list(iopbend)%a + first_atom - 1
507 25 : index_b = opbend_list(iopbend)%b + first_atom - 1
508 25 : index_c = opbend_list(iopbend)%c + first_atom - 1
509 25 : index_d = opbend_list(iopbend)%d + first_atom - 1
510 :
511 100 : t12 = particle_set(index_a)%r - particle_set(index_b)%r
512 100 : t32 = particle_set(index_c)%r - particle_set(index_b)%r
513 100 : t34 = particle_set(index_c)%r - particle_set(index_d)%r
514 100 : t43 = particle_set(index_d)%r - particle_set(index_c)%r
515 100 : t41 = particle_set(index_d)%r - particle_set(index_a)%r
516 100 : t42 = pbc(t41 + t12, cell)
517 100 : t12 = pbc(t12, cell)
518 100 : t32 = pbc(t32, cell)
519 100 : t41 = pbc(t41, cell)
520 100 : t43 = pbc(t43, cell)
521 : ! tm = t32 x t12
522 25 : tm(1) = t32(2)*t12(3) - t12(2)*t32(3)
523 25 : tm(2) = -t32(1)*t12(3) + t12(1)*t32(3)
524 25 : tm(3) = t32(1)*t12(2) - t12(1)*t32(2)
525 100 : sm = SQRT(DOT_PRODUCT(tm, tm))
526 100 : 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 25 : gt1, gt2, gt3, gt4, energy, fscalar)
532 25 : pot_opbend = pot_opbend + energy
533 25 : IF (atener) THEN
534 0 : ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
535 0 : ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
536 0 : ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
537 0 : ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
538 : END IF
539 25 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
540 25 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
541 25 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
542 25 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
543 25 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
544 25 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
545 25 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
546 25 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
547 25 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
548 25 : particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
549 25 : particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
550 25 : particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
551 :
552 : ! computing the pressure tensor
553 100 : k1 = fscalar*gt1
554 100 : k3 = fscalar*gt3
555 100 : k4 = fscalar*gt4
556 :
557 25 : IF (use_virial) CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_opbend)
558 25 : IF (atstress) THEN
559 0 : k1 = 0.25_dp*fscalar*gt1
560 0 : k3 = 0.25_dp*fscalar*gt3
561 0 : k4 = 0.25_dp*fscalar*gt4
562 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_a))
563 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_b))
564 0 : CALL get_pv_torsion(k1, k3, k4, t12, t32, t43, pv_a(:, :, index_c))
565 0 : 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 2050618 : IF (PRESENT(f_opbend)) THEN
570 0 : f_opbend(:, index_a) = f_opbend(:, index_a) + fscalar*gt1
571 0 : f_opbend(:, index_b) = f_opbend(:, index_b) + fscalar*gt2
572 0 : f_opbend(:, index_c) = f_opbend(:, index_c) + fscalar*gt3
573 0 : 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 83088 : CALL timestop(handle)
580 :
581 83088 : END SUBROUTINE force_intra_control
582 :
583 : END MODULE fist_intra_force
584 :
|
