Line data Source code
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 : ! **************************************************************************************************
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 158304 : 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 79152 : pot_opbend, pot_shell, pv_bond, pv_bend, pv_urey_bradley, pv_torsion, &
85 237456 : pv_imp_torsion, pv_opbend, f_bond, f_bend, f_torsion, f_ub, &
86 158304 : 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 79152 : REAL(KIND=dp), DIMENSION(:), POINTER :: ener_a
119 79152 : TYPE(bend_type), POINTER :: bend_list(:)
120 79152 : TYPE(bond_type), POINTER :: bond_list(:)
121 : TYPE(cp_logger_type), POINTER :: logger
122 79152 : TYPE(impr_type), POINTER :: impr_list(:)
123 : TYPE(molecule_kind_type), POINTER :: molecule_kind
124 : TYPE(molecule_type), POINTER :: molecule
125 79152 : TYPE(opbend_type), POINTER :: opbend_list(:)
126 79152 : TYPE(shell_type), POINTER :: shell_list(:)
127 79152 : TYPE(torsion_type), POINTER :: torsion_list(:)
128 79152 : TYPE(ub_type), POINTER :: ub_list(:)
129 :
130 79152 : CALL timeset(routineN, handle)
131 79152 : NULLIFY (logger)
132 79152 : logger => cp_get_default_logger()
133 :
134 79152 : IF (PRESENT(f_bond)) f_bond = 0.0_dp
135 79152 : IF (PRESENT(f_bend)) f_bend = 0.0_dp
136 79152 : IF (PRESENT(f_torsion)) f_torsion = 0.0_dp
137 79152 : IF (PRESENT(f_imptor)) f_imptor = 0.0_dp
138 79152 : IF (PRESENT(f_ub)) f_ub = 0.0_dp
139 :
140 79152 : pot_bond = 0.0_dp
141 79152 : pot_bend = 0.0_dp
142 79152 : pot_urey_bradley = 0.0_dp
143 79152 : pot_torsion = 0.0_dp
144 79152 : pot_imp_torsion = 0.0_dp
145 79152 : pot_opbend = 0.0_dp
146 79152 : pot_shell = 0.0_dp
147 :
148 79152 : atener = atprop_env%energy
149 79152 : IF (atener) ener_a => atprop_env%atener
150 :
151 79152 : nkind = SIZE(molecule_kind_set)
152 1402490 : MOL: DO ikind = 1, nkind
153 1323338 : nmol_per_kind = local_molecules%n_el(ikind)
154 :
155 3431372 : DO imol = 1, nmol_per_kind
156 2028882 : i = local_molecules%list(ikind)%array(imol)
157 2028882 : molecule => molecule_set(i)
158 2028882 : 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 2028882 : shell_list=shell_list)
166 :
167 2028882 : CALL get_molecule(molecule, first_atom=first_atom)
168 :
169 4802673 : BOND: DO ibond = 1, nbonds
170 2773791 : index_a = bond_list(ibond)%a + first_atom - 1
171 2773791 : index_b = bond_list(ibond)%b + first_atom - 1
172 11095164 : rij = particle_set(index_a)%r - particle_set(index_b)%r
173 11095164 : 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 2773791 : energy, fscalar)
179 2773791 : pot_bond = pot_bond + energy
180 2773791 : IF (atener) THEN
181 606 : ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
182 606 : ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
183 : END IF
184 :
185 2773791 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
186 2773791 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
187 2773791 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
188 2773791 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
189 2773791 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
190 2773791 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
191 :
192 : ! computing the pressure tensor
193 11095164 : k2 = -rij*fscalar
194 2773791 : IF (use_virial) CALL get_pv_bond(k2, rij, pv_bond)
195 :
196 : ! the contribution from the bonds. ONLY FOR DEBUG
197 4802673 : IF (PRESENT(f_bond)) THEN
198 0 : f_bond(1, index_a) = f_bond(1, index_a) - rij(1)*fscalar
199 0 : f_bond(2, index_a) = f_bond(2, index_a) - rij(2)*fscalar
200 0 : f_bond(3, index_a) = f_bond(3, index_a) - rij(3)*fscalar
201 0 : f_bond(1, index_b) = f_bond(1, index_b) + rij(1)*fscalar
202 0 : f_bond(2, index_b) = f_bond(2, index_b) + rij(2)*fscalar
203 0 : f_bond(3, index_b) = f_bond(3, index_b) + rij(3)*fscalar
204 : END IF
205 :
206 : END DO BOND
207 :
208 2445279 : SHELL: DO ishell = 1, nshell
209 416397 : index_a = shell_list(ishell)%a + first_atom - 1
210 416397 : index_b = particle_set(index_a)%shell_index
211 1665588 : rij = core_particle_set(index_b)%r - shell_particle_set(index_b)%r
212 1665588 : rij = pbc(rij, cell)
213 416397 : k2_spring = shell_list(ishell)%shell_kind%k2_spring
214 416397 : k4_spring = shell_list(ishell)%shell_kind%k4_spring
215 1665588 : r2 = DOT_PRODUCT(rij, rij)
216 416397 : energy = 0.5_dp*(k2_spring + k4_spring*r2/12.0_dp)*r2
217 416397 : fscalar = k2_spring + k4_spring*r2/6.0_dp
218 416397 : pot_shell = pot_shell + energy
219 416397 : IF (atener) THEN
220 1080 : ener_a(index_a) = ener_a(index_a) + energy
221 : END IF
222 416397 : core_particle_set(index_b)%f(1) = core_particle_set(index_b)%f(1) - rij(1)*fscalar
223 416397 : core_particle_set(index_b)%f(2) = core_particle_set(index_b)%f(2) - rij(2)*fscalar
224 416397 : core_particle_set(index_b)%f(3) = core_particle_set(index_b)%f(3) - rij(3)*fscalar
225 416397 : shell_particle_set(index_b)%f(1) = shell_particle_set(index_b)%f(1) + rij(1)*fscalar
226 416397 : shell_particle_set(index_b)%f(2) = shell_particle_set(index_b)%f(2) + rij(2)*fscalar
227 416397 : 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 2445279 : IF (use_virial) THEN
230 1369328 : k1 = -rij*fscalar
231 342332 : CALL get_pv_bond(k1, rij, pv_bond)
232 : END IF
233 : END DO SHELL
234 :
235 2137085 : UREY_BRADLEY: DO ibend = 1, nub
236 108203 : index_a = ub_list(ibend)%a + first_atom - 1
237 108203 : index_b = ub_list(ibend)%c + first_atom - 1
238 432812 : rij = particle_set(index_a)%r - particle_set(index_b)%r
239 432812 : rij = pbc(rij, cell)
240 : CALL force_bonds(ub_list(ibend)%ub_kind%id_type, rij, &
241 : ub_list(ibend)%ub_kind%r0, &
242 108203 : ub_list(ibend)%ub_kind%k, 0.0_dp, energy, fscalar)
243 108203 : pot_urey_bradley = pot_urey_bradley + energy
244 108203 : IF (atener) THEN
245 0 : ener_a(index_a) = ener_a(index_a) + 0.5_dp*energy
246 0 : ener_a(index_b) = ener_a(index_b) + 0.5_dp*energy
247 : END IF
248 108203 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) - rij(1)*fscalar
249 108203 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) - rij(2)*fscalar
250 108203 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) - rij(3)*fscalar
251 108203 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + rij(1)*fscalar
252 108203 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + rij(2)*fscalar
253 108203 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + rij(3)*fscalar
254 :
255 : ! computing the pressure tensor
256 432812 : k2 = -rij*fscalar
257 108203 : IF (use_virial) CALL get_pv_bond(k2, rij, pv_urey_bradley)
258 :
259 : ! the contribution from the ub. ONLY FOR DEBUG
260 2137085 : IF (PRESENT(f_ub)) THEN
261 0 : f_ub(:, index_a) = f_ub(:, index_a) - rij*fscalar
262 0 : f_ub(:, index_b) = f_ub(:, index_b) + rij*fscalar
263 : END IF
264 :
265 : END DO UREY_BRADLEY
266 :
267 3380462 : BEND: DO ibend = 1, nbends
268 1351580 : index_a = bend_list(ibend)%a + first_atom - 1
269 1351580 : index_b = bend_list(ibend)%b + first_atom - 1
270 1351580 : index_c = bend_list(ibend)%c + first_atom - 1
271 5406320 : b12 = particle_set(index_a)%r - particle_set(index_b)%r
272 5406320 : b32 = particle_set(index_c)%r - particle_set(index_b)%r
273 5406320 : b12 = pbc(b12, cell)
274 5406320 : b32 = pbc(b32, cell)
275 5406320 : d12 = SQRT(DOT_PRODUCT(b12, b12))
276 1351580 : id12 = 1.0_dp/d12
277 5406320 : d32 = SQRT(DOT_PRODUCT(b32, b32))
278 1351580 : id32 = 1.0_dp/d32
279 5406320 : dist = DOT_PRODUCT(b12, b32)
280 1351580 : theta = (dist*id12*id32)
281 1351580 : IF (theta < -1.0_dp) theta = -1.0_dp
282 1351580 : IF (theta > +1.0_dp) theta = +1.0_dp
283 1351580 : 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 1351580 : g1, g2, g3, energy, fscalar)
296 1351580 : pot_bend = pot_bend + energy
297 1351580 : IF (atener) THEN
298 303 : ener_a(index_a) = ener_a(index_a) + energy/3._dp
299 303 : ener_a(index_b) = ener_a(index_b) + energy/3._dp
300 303 : ener_a(index_c) = ener_a(index_c) + energy/3._dp
301 : END IF
302 1351580 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + g1(1)*fscalar
303 1351580 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + g1(2)*fscalar
304 1351580 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + g1(3)*fscalar
305 1351580 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + g2(1)*fscalar
306 1351580 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + g2(2)*fscalar
307 1351580 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + g2(3)*fscalar
308 1351580 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + g3(1)*fscalar
309 1351580 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + g3(2)*fscalar
310 1351580 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + g3(3)*fscalar
311 :
312 : ! computing the pressure tensor
313 5406320 : k1 = fscalar*g1
314 5406320 : k3 = fscalar*g3
315 1351580 : IF (use_virial) CALL get_pv_bend(k1, k3, b12, b32, pv_bend)
316 :
317 : ! the contribution from the bends. ONLY FOR DEBUG
318 3380462 : IF (PRESENT(f_bend)) THEN
319 0 : f_bend(:, index_a) = f_bend(:, index_a) + fscalar*g1
320 0 : f_bend(:, index_b) = f_bend(:, index_b) + fscalar*g2
321 0 : f_bend(:, index_c) = f_bend(:, index_c) + fscalar*g3
322 : END IF
323 : END DO BEND
324 :
325 2716188 : TORSION: DO itorsion = 1, ntorsions
326 687306 : index_a = torsion_list(itorsion)%a + first_atom - 1
327 687306 : index_b = torsion_list(itorsion)%b + first_atom - 1
328 687306 : index_c = torsion_list(itorsion)%c + first_atom - 1
329 687306 : index_d = torsion_list(itorsion)%d + first_atom - 1
330 2749224 : t12 = particle_set(index_a)%r - particle_set(index_b)%r
331 2749224 : t32 = particle_set(index_c)%r - particle_set(index_b)%r
332 2749224 : t34 = particle_set(index_c)%r - particle_set(index_d)%r
333 2749224 : t43 = particle_set(index_d)%r - particle_set(index_c)%r
334 2749224 : t12 = pbc(t12, cell)
335 2749224 : t32 = pbc(t32, cell)
336 2749224 : t34 = pbc(t34, cell)
337 2749224 : t43 = pbc(t43, cell)
338 : ! t12 x t32
339 687306 : tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
340 687306 : tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
341 687306 : tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
342 : ! t32 x t34
343 687306 : tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
344 687306 : tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
345 687306 : tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
346 2749224 : sm = SQRT(DOT_PRODUCT(tm, tm))
347 687306 : ism = 1.0_dp/sm
348 2749224 : sn = SQRT(DOT_PRODUCT(tn, tn))
349 687306 : isn = 1.0_dp/sn
350 2749224 : s32 = SQRT(DOT_PRODUCT(t32, t32))
351 687306 : is32 = 1.0_dp/s32
352 2749224 : dist1 = DOT_PRODUCT(t12, t32)
353 2749224 : dist2 = DOT_PRODUCT(t34, t32)
354 3462211 : 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 746023 : gt1, gt2, gt3, gt4, energy, fscalar)
361 746023 : pot_torsion = pot_torsion + energy
362 746023 : IF (atener) THEN
363 0 : ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
364 0 : ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
365 0 : ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
366 0 : ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
367 : END IF
368 746023 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
369 746023 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
370 746023 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
371 746023 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
372 746023 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
373 746023 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
374 746023 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
375 746023 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
376 746023 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
377 746023 : particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
378 746023 : particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
379 746023 : particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
380 :
381 : ! computing the pressure tensor
382 2984092 : k1 = fscalar*gt1
383 2984092 : k3 = fscalar*gt3
384 2984092 : k4 = fscalar*gt4
385 746023 : 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 1433329 : IF (PRESENT(f_torsion)) THEN
389 0 : f_torsion(:, index_a) = f_torsion(:, index_a) + fscalar*gt1
390 0 : f_torsion(:, index_b) = f_torsion(:, index_b) + fscalar*gt2
391 0 : f_torsion(:, index_c) = f_torsion(:, index_c) + fscalar*gt3
392 0 : f_torsion(:, index_d) = f_torsion(:, index_d) + fscalar*gt4
393 : END IF
394 : END DO
395 : END DO TORSION
396 :
397 2048189 : IMP_TORSION: DO itorsion = 1, nimptors
398 19307 : index_a = impr_list(itorsion)%a + first_atom - 1
399 19307 : index_b = impr_list(itorsion)%b + first_atom - 1
400 19307 : index_c = impr_list(itorsion)%c + first_atom - 1
401 19307 : index_d = impr_list(itorsion)%d + first_atom - 1
402 77228 : t12 = particle_set(index_a)%r - particle_set(index_b)%r
403 77228 : t32 = particle_set(index_c)%r - particle_set(index_b)%r
404 77228 : t34 = particle_set(index_c)%r - particle_set(index_d)%r
405 77228 : t43 = particle_set(index_d)%r - particle_set(index_c)%r
406 77228 : t12 = pbc(t12, cell)
407 77228 : t32 = pbc(t32, cell)
408 77228 : t34 = pbc(t34, cell)
409 77228 : t43 = pbc(t43, cell)
410 : ! t12 x t32
411 19307 : tm(1) = t12(2)*t32(3) - t32(2)*t12(3)
412 19307 : tm(2) = -t12(1)*t32(3) + t32(1)*t12(3)
413 19307 : tm(3) = t12(1)*t32(2) - t32(1)*t12(2)
414 : ! t32 x t34
415 19307 : tn(1) = t32(2)*t34(3) - t34(2)*t32(3)
416 19307 : tn(2) = -t32(1)*t34(3) + t34(1)*t32(3)
417 19307 : tn(3) = t32(1)*t34(2) - t34(1)*t32(2)
418 77228 : sm = SQRT(DOT_PRODUCT(tm, tm))
419 19307 : ism = 1.0_dp/sm
420 77228 : sn = SQRT(DOT_PRODUCT(tn, tn))
421 19307 : isn = 1.0_dp/sn
422 77228 : s32 = SQRT(DOT_PRODUCT(t32, t32))
423 19307 : is32 = 1.0_dp/s32
424 77228 : dist1 = DOT_PRODUCT(t12, t32)
425 77228 : 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 19307 : gt1, gt2, gt3, gt4, energy, fscalar)
431 19307 : pot_imp_torsion = pot_imp_torsion + energy
432 19307 : IF (atener) THEN
433 0 : ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
434 0 : ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
435 0 : ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
436 0 : ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
437 : END IF
438 19307 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
439 19307 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
440 19307 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
441 19307 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
442 19307 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
443 19307 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
444 19307 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
445 19307 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
446 19307 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
447 19307 : particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
448 19307 : particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
449 19307 : particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
450 :
451 : ! computing the pressure tensor
452 77228 : k1 = fscalar*gt1
453 77228 : k3 = fscalar*gt3
454 77228 : k4 = fscalar*gt4
455 19307 : 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 2048189 : IF (PRESENT(f_imptor)) THEN
459 0 : f_imptor(:, index_a) = f_imptor(:, index_a) + fscalar*gt1
460 0 : f_imptor(:, index_b) = f_imptor(:, index_b) + fscalar*gt2
461 0 : f_imptor(:, index_c) = f_imptor(:, index_c) + fscalar*gt3
462 0 : f_imptor(:, index_d) = f_imptor(:, index_d) + fscalar*gt4
463 : END IF
464 : END DO IMP_TORSION
465 :
466 5381127 : OPBEND: DO iopbend = 1, nopbends
467 25 : index_a = opbend_list(iopbend)%a + first_atom - 1
468 25 : index_b = opbend_list(iopbend)%b + first_atom - 1
469 25 : index_c = opbend_list(iopbend)%c + first_atom - 1
470 25 : index_d = opbend_list(iopbend)%d + first_atom - 1
471 :
472 100 : t12 = particle_set(index_a)%r - particle_set(index_b)%r
473 100 : t32 = particle_set(index_c)%r - particle_set(index_b)%r
474 100 : t34 = particle_set(index_c)%r - particle_set(index_d)%r
475 100 : t43 = particle_set(index_d)%r - particle_set(index_c)%r
476 100 : t41 = particle_set(index_d)%r - particle_set(index_a)%r
477 100 : t42 = pbc(t41 + t12, cell)
478 100 : t12 = pbc(t12, cell)
479 100 : t32 = pbc(t32, cell)
480 100 : t41 = pbc(t41, cell)
481 100 : t43 = pbc(t43, cell)
482 : ! tm = t32 x t12
483 25 : tm(1) = t32(2)*t12(3) - t12(2)*t32(3)
484 25 : tm(2) = -t32(1)*t12(3) + t12(1)*t32(3)
485 25 : tm(3) = t32(1)*t12(2) - t12(1)*t32(2)
486 100 : sm = SQRT(DOT_PRODUCT(tm, tm))
487 100 : 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 25 : gt1, gt2, gt3, gt4, energy, fscalar)
493 25 : pot_opbend = pot_opbend + energy
494 25 : IF (atener) THEN
495 0 : ener_a(index_a) = ener_a(index_a) + energy*0.25_dp
496 0 : ener_a(index_b) = ener_a(index_b) + energy*0.25_dp
497 0 : ener_a(index_c) = ener_a(index_c) + energy*0.25_dp
498 0 : ener_a(index_d) = ener_a(index_d) + energy*0.25_dp
499 : END IF
500 25 : particle_set(index_a)%f(1) = particle_set(index_a)%f(1) + gt1(1)*fscalar
501 25 : particle_set(index_a)%f(2) = particle_set(index_a)%f(2) + gt1(2)*fscalar
502 25 : particle_set(index_a)%f(3) = particle_set(index_a)%f(3) + gt1(3)*fscalar
503 25 : particle_set(index_b)%f(1) = particle_set(index_b)%f(1) + gt2(1)*fscalar
504 25 : particle_set(index_b)%f(2) = particle_set(index_b)%f(2) + gt2(2)*fscalar
505 25 : particle_set(index_b)%f(3) = particle_set(index_b)%f(3) + gt2(3)*fscalar
506 25 : particle_set(index_c)%f(1) = particle_set(index_c)%f(1) + gt3(1)*fscalar
507 25 : particle_set(index_c)%f(2) = particle_set(index_c)%f(2) + gt3(2)*fscalar
508 25 : particle_set(index_c)%f(3) = particle_set(index_c)%f(3) + gt3(3)*fscalar
509 25 : particle_set(index_d)%f(1) = particle_set(index_d)%f(1) + gt4(1)*fscalar
510 25 : particle_set(index_d)%f(2) = particle_set(index_d)%f(2) + gt4(2)*fscalar
511 25 : particle_set(index_d)%f(3) = particle_set(index_d)%f(3) + gt4(3)*fscalar
512 :
513 : ! computing the pressure tensor
514 100 : k1 = fscalar*gt1
515 100 : k3 = fscalar*gt3
516 100 : k4 = fscalar*gt4
517 :
518 25 : 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 2028907 : IF (PRESENT(f_opbend)) THEN
522 0 : f_opbend(:, index_a) = f_opbend(:, index_a) + fscalar*gt1
523 0 : f_opbend(:, index_b) = f_opbend(:, index_b) + fscalar*gt2
524 0 : f_opbend(:, index_c) = f_opbend(:, index_c) + fscalar*gt3
525 0 : 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 79152 : CALL timestop(handle)
532 :
533 79152 : END SUBROUTINE force_intra_control
534 :
535 : END MODULE fist_intra_force
536 :
|
