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 : !> Splitting and cleaning the original force_field_pack - May 2007
11 : !> Teodoro Laino - Zurich University
12 : !> \author CJM
13 : ! **************************************************************************************************
14 : MODULE force_fields_all
15 :
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind,&
18 : get_atomic_kind_set,&
19 : set_atomic_kind
20 : USE atoms_input, ONLY: read_shell_coord_input
21 : USE cell_types, ONLY: cell_type
22 : USE cp_linked_list_input, ONLY: cp_sll_val_next,&
23 : cp_sll_val_type
24 : USE cp_log_handling, ONLY: cp_to_string
25 : USE damping_dipole_types, ONLY: damping_p_create,&
26 : damping_p_type,&
27 : tang_toennies
28 : USE ewald_environment_types, ONLY: ewald_env_get,&
29 : ewald_env_set,&
30 : ewald_environment_type
31 : USE external_potential_types, ONLY: fist_potential_type,&
32 : get_potential,&
33 : set_potential
34 : USE force_field_kind_types, ONLY: &
35 : allocate_bend_kind_set, allocate_bond_kind_set, allocate_impr_kind_set, &
36 : allocate_opbend_kind_set, allocate_torsion_kind_set, allocate_ub_kind_set, bend_kind_type, &
37 : bond_kind_type, do_ff_amber, do_ff_charmm, do_ff_g87, do_ff_g96, do_ff_undef, &
38 : impr_kind_type, opbend_kind_type, torsion_kind_type, ub_kind_type
39 : USE force_field_types, ONLY: amber_info_type,&
40 : charmm_info_type,&
41 : force_field_type,&
42 : gromos_info_type,&
43 : input_info_type
44 : USE input_constants, ONLY: do_qmmm_none
45 : USE input_cp2k_binary_restarts, ONLY: read_binary_cs_coordinates
46 : USE input_section_types, ONLY: section_vals_get,&
47 : section_vals_get_subs_vals,&
48 : section_vals_list_get,&
49 : section_vals_type,&
50 : section_vals_val_get
51 : USE input_val_types, ONLY: val_get,&
52 : val_type
53 : USE kinds, ONLY: default_path_length,&
54 : default_string_length,&
55 : dp
56 : USE mathconstants, ONLY: sqrthalf
57 : USE memory_utilities, ONLY: reallocate
58 : USE molecule_kind_types, ONLY: &
59 : bend_type, bond_type, get_molecule_kind, impr_type, molecule_kind_type, opbend_type, &
60 : set_molecule_kind, shell_type, torsion_type, ub_type
61 : USE molecule_types, ONLY: get_molecule,&
62 : molecule_type
63 : USE pair_potential, ONLY: get_nonbond_storage,&
64 : spline_nonbond_control
65 : USE pair_potential_coulomb, ONLY: potential_coulomb
66 : USE pair_potential_types, ONLY: &
67 : ace_type, allegro_type, deepmd_type, ea_type, lj_charmm_type, lj_type, nequip_type, &
68 : nn_type, nosh_nosh, nosh_sh, pair_potential_lj_create, pair_potential_pp_create, &
69 : pair_potential_pp_type, pair_potential_single_add, pair_potential_single_clean, &
70 : pair_potential_single_copy, pair_potential_single_type, quip_type, sh_sh, siepmann_type, &
71 : tersoff_type
72 : USE particle_types, ONLY: allocate_particle_set,&
73 : particle_type
74 : USE physcon, ONLY: bohr
75 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
76 : USE qmmm_types_low, ONLY: qmmm_env_mm_type
77 : USE shell_potential_types, ONLY: shell_kind_type
78 : USE splines_types, ONLY: spline_data_p_release,&
79 : spline_data_p_retain,&
80 : spline_data_p_type,&
81 : spline_env_release,&
82 : spline_environment_type
83 : USE string_utilities, ONLY: compress,&
84 : integer_to_string,&
85 : uppercase
86 : #include "./base/base_uses.f90"
87 :
88 : IMPLICIT NONE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_all'
91 :
92 : PRIVATE
93 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
94 :
95 : PUBLIC :: force_field_unique_bond, &
96 : force_field_unique_bend, &
97 : force_field_unique_ub, &
98 : force_field_unique_tors, &
99 : force_field_unique_impr, &
100 : force_field_unique_opbend, &
101 : force_field_pack_bond, &
102 : force_field_pack_bend, &
103 : force_field_pack_ub, &
104 : force_field_pack_tors, &
105 : force_field_pack_impr, &
106 : force_field_pack_opbend, &
107 : force_field_pack_charge, &
108 : force_field_pack_charges, &
109 : force_field_pack_radius, &
110 : force_field_pack_pol, &
111 : force_field_pack_shell, &
112 : force_field_pack_nonbond14, &
113 : force_field_pack_nonbond, &
114 : force_field_pack_splines, &
115 : force_field_pack_eicut, &
116 : force_field_pack_damp
117 :
118 : CONTAINS
119 :
120 : ! **************************************************************************************************
121 : !> \brief Determine the number of unique bond kind and allocate bond_kind_set
122 : !> \param particle_set ...
123 : !> \param molecule_kind_set ...
124 : !> \param molecule_set ...
125 : !> \param ff_type ...
126 : !> \param iw ...
127 : ! **************************************************************************************************
128 2643 : SUBROUTINE force_field_unique_bond(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
129 :
130 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
131 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
132 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
133 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
134 : INTEGER, INTENT(IN) :: iw
135 :
136 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bond'
137 :
138 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, &
139 : name_atm_b2
140 : INTEGER :: atm_a, atm_b, counter, first, handle2, &
141 : i, j, k, last, natom, nbond
142 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
143 2643 : INTEGER, POINTER :: map_bond_kind(:)
144 : LOGICAL :: found
145 : TYPE(atomic_kind_type), POINTER :: atomic_kind
146 2643 : TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set
147 2643 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
148 : TYPE(molecule_kind_type), POINTER :: molecule_kind
149 : TYPE(molecule_type), POINTER :: molecule
150 :
151 2643 : CALL timeset(routineN, handle2)
152 :
153 2643 : IF (iw > 0) THEN
154 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
155 239 : "FORCEFIELD| Checking for unique bond terms"
156 : END IF
157 :
158 75465 : DO i = 1, SIZE(molecule_kind_set)
159 72822 : molecule_kind => molecule_kind_set(i)
160 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
161 : molecule_list=molecule_list, &
162 : natom=natom, &
163 72822 : nbond=nbond, bond_list=bond_list)
164 72822 : molecule => molecule_set(molecule_list(1))
165 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
166 148287 : IF (nbond > 0) THEN
167 88287 : ALLOCATE (map_bond_kind(nbond))
168 29429 : counter = 0
169 29429 : IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN
170 148 : DO j = 1, nbond
171 148 : map_bond_kind(j) = j
172 : END DO
173 20 : counter = nbond
174 : ELSE
175 143456 : DO j = 1, nbond
176 114047 : atm_a = bond_list(j)%a
177 114047 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
178 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
179 114047 : name=name_atm_a)
180 114047 : atm_b = bond_list(j)%b
181 114047 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
182 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
183 114047 : name=name_atm_b)
184 114047 : found = .FALSE.
185 482677 : DO k = 1, j - 1
186 415898 : atm_a = bond_list(k)%a
187 415898 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
188 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
189 415898 : name=name_atm_a2)
190 415898 : atm_b = bond_list(k)%b
191 415898 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
192 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
193 415898 : name=name_atm_b2)
194 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
195 415898 : ((name_atm_b) == (name_atm_b2))) .OR. &
196 : (((name_atm_a) == (name_atm_b2)) .AND. &
197 66779 : ((name_atm_b) == (name_atm_a2)))) THEN
198 47268 : found = .TRUE.
199 47268 : map_bond_kind(j) = map_bond_kind(k)
200 : EXIT
201 : END IF
202 : END DO
203 29409 : IF (.NOT. found) THEN
204 66779 : counter = counter + 1
205 66779 : map_bond_kind(j) = counter
206 : END IF
207 : END DO
208 : END IF
209 29429 : NULLIFY (bond_kind_set)
210 29429 : CALL allocate_bond_kind_set(bond_kind_set, counter)
211 143604 : DO j = 1, nbond
212 143604 : bond_list(j)%bond_kind => bond_kind_set(map_bond_kind(j))
213 : END DO
214 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
215 29429 : bond_kind_set=bond_kind_set, bond_list=bond_list)
216 29429 : DEALLOCATE (map_bond_kind)
217 : END IF
218 : END DO
219 2643 : CALL timestop(handle2)
220 :
221 2643 : END SUBROUTINE force_field_unique_bond
222 :
223 : ! **************************************************************************************************
224 : !> \brief Determine the number of unique bend kind and allocate bend_kind_set
225 : !> \param particle_set ...
226 : !> \param molecule_kind_set ...
227 : !> \param molecule_set ...
228 : !> \param ff_type ...
229 : !> \param iw ...
230 : ! **************************************************************************************************
231 2643 : SUBROUTINE force_field_unique_bend(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
232 :
233 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
234 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
235 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
236 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
237 : INTEGER, INTENT(IN) :: iw
238 :
239 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bend'
240 :
241 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, &
242 : name_atm_b2, name_atm_c, name_atm_c2
243 : INTEGER :: atm_a, atm_b, atm_c, counter, first, &
244 : handle2, i, j, k, last, natom, nbend
245 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
246 2643 : INTEGER, POINTER :: map_bend_kind(:)
247 : LOGICAL :: found
248 : TYPE(atomic_kind_type), POINTER :: atomic_kind
249 2643 : TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set
250 2643 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
251 : TYPE(molecule_kind_type), POINTER :: molecule_kind
252 : TYPE(molecule_type), POINTER :: molecule
253 :
254 2643 : CALL timeset(routineN, handle2)
255 :
256 2643 : IF (iw > 0) THEN
257 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
258 239 : "FORCEFIELD| Checking for unique bend terms"
259 : END IF
260 :
261 75465 : DO i = 1, SIZE(molecule_kind_set)
262 72822 : molecule_kind => molecule_kind_set(i)
263 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
264 : molecule_list=molecule_list, &
265 : natom=natom, &
266 72822 : nbend=nbend, bend_list=bend_list)
267 72822 : molecule => molecule_set(molecule_list(1))
268 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
269 148287 : IF (nbend > 0) THEN
270 87297 : ALLOCATE (map_bend_kind(nbend))
271 29099 : counter = 0
272 29099 : IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN
273 168 : DO j = 1, nbend
274 168 : map_bend_kind(j) = j
275 : END DO
276 12 : counter = nbend
277 : ELSE
278 168025 : DO j = 1, nbend
279 138938 : atm_a = bend_list(j)%a
280 138938 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
281 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
282 138938 : name=name_atm_a)
283 138938 : atm_b = bend_list(j)%b
284 138938 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
285 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
286 138938 : name=name_atm_b)
287 138938 : atm_c = bend_list(j)%c
288 138938 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
289 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
290 138938 : name=name_atm_c)
291 138938 : found = .FALSE.
292 2277017 : DO k = 1, j - 1
293 2182221 : atm_a = bend_list(k)%a
294 2182221 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
295 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
296 2182221 : name=name_atm_a2)
297 2182221 : atm_b = bend_list(k)%b
298 2182221 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
299 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
300 2182221 : name=name_atm_b2)
301 2182221 : atm_c = bend_list(k)%c
302 2182221 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
303 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
304 2182221 : name=name_atm_c2)
305 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
306 : ((name_atm_b) == (name_atm_b2)) .AND. &
307 2182221 : ((name_atm_c) == (name_atm_c2))) .OR. &
308 : (((name_atm_a) == (name_atm_c2)) .AND. &
309 : ((name_atm_b) == (name_atm_b2)) .AND. &
310 94796 : ((name_atm_c) == (name_atm_a2)))) THEN
311 44142 : found = .TRUE.
312 44142 : map_bend_kind(j) = map_bend_kind(k)
313 : EXIT
314 : END IF
315 : END DO
316 29087 : IF (.NOT. found) THEN
317 94796 : counter = counter + 1
318 94796 : map_bend_kind(j) = counter
319 : END IF
320 : END DO
321 : END IF
322 29099 : NULLIFY (bend_kind_set)
323 29099 : CALL allocate_bend_kind_set(bend_kind_set, counter)
324 168193 : DO j = 1, nbend
325 168193 : bend_list(j)%bend_kind => bend_kind_set(map_bend_kind(j))
326 : END DO
327 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
328 29099 : bend_kind_set=bend_kind_set, bend_list=bend_list)
329 29099 : DEALLOCATE (map_bend_kind)
330 : END IF
331 : END DO
332 :
333 2643 : CALL timestop(handle2)
334 :
335 2643 : END SUBROUTINE force_field_unique_bend
336 :
337 : ! **************************************************************************************************
338 : !> \brief Determine the number of unique Urey-Bradley kind and allocate ub_kind_set
339 : !> \param particle_set ...
340 : !> \param molecule_kind_set ...
341 : !> \param molecule_set ...
342 : !> \param iw ...
343 : ! **************************************************************************************************
344 2643 : SUBROUTINE force_field_unique_ub(particle_set, molecule_kind_set, molecule_set, iw)
345 :
346 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
347 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
348 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
349 : INTEGER, INTENT(IN) :: iw
350 :
351 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_ub'
352 :
353 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, &
354 : name_atm_b2, name_atm_c, name_atm_c2
355 : INTEGER :: atm_a, atm_b, atm_c, counter, first, &
356 : handle2, i, j, k, last, natom, nub
357 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
358 2643 : INTEGER, POINTER :: map_ub_kind(:)
359 : LOGICAL :: found
360 : TYPE(atomic_kind_type), POINTER :: atomic_kind
361 : TYPE(molecule_kind_type), POINTER :: molecule_kind
362 : TYPE(molecule_type), POINTER :: molecule
363 2643 : TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set
364 2643 : TYPE(ub_type), DIMENSION(:), POINTER :: ub_list
365 :
366 2643 : CALL timeset(routineN, handle2)
367 :
368 2643 : IF (iw > 0) THEN
369 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
370 239 : "FORCEFIELD| Checking for unique Urey-Bradley terms"
371 : END IF
372 :
373 75465 : DO i = 1, SIZE(molecule_kind_set)
374 72822 : molecule_kind => molecule_kind_set(i)
375 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
376 : molecule_list=molecule_list, &
377 : natom=natom, &
378 72822 : nub=nub, ub_list=ub_list)
379 72822 : molecule => molecule_set(molecule_list(1))
380 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
381 148287 : IF (nub > 0) THEN
382 87255 : ALLOCATE (map_ub_kind(nub))
383 29085 : counter = 0
384 168021 : DO j = 1, nub
385 138936 : atm_a = ub_list(j)%a
386 138936 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
387 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
388 138936 : name=name_atm_a)
389 138936 : atm_b = ub_list(j)%b
390 138936 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
391 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
392 138936 : name=name_atm_b)
393 138936 : atm_c = ub_list(j)%c
394 138936 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
395 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
396 138936 : name=name_atm_c)
397 138936 : found = .FALSE.
398 2277015 : DO k = 1, j - 1
399 2182221 : atm_a = ub_list(k)%a
400 2182221 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
401 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
402 2182221 : name=name_atm_a2)
403 2182221 : atm_b = ub_list(k)%b
404 2182221 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
405 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
406 2182221 : name=name_atm_b2)
407 2182221 : atm_c = ub_list(k)%c
408 2182221 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
409 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
410 2182221 : name=name_atm_c2)
411 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
412 : ((name_atm_b) == (name_atm_b2)) .AND. &
413 2182221 : ((name_atm_c) == (name_atm_c2))) .OR. &
414 : (((name_atm_a) == (name_atm_c2)) .AND. &
415 : ((name_atm_b) == (name_atm_b2)) .AND. &
416 94794 : ((name_atm_c) == (name_atm_a2)))) THEN
417 44142 : found = .TRUE.
418 44142 : map_ub_kind(j) = map_ub_kind(k)
419 : EXIT
420 : END IF
421 : END DO
422 29085 : IF (.NOT. found) THEN
423 94794 : counter = counter + 1
424 94794 : map_ub_kind(j) = counter
425 : END IF
426 : END DO
427 29085 : CALL allocate_ub_kind_set(ub_kind_set, counter)
428 168021 : DO j = 1, nub
429 168021 : ub_list(j)%ub_kind => ub_kind_set(map_ub_kind(j))
430 : END DO
431 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
432 29085 : ub_kind_set=ub_kind_set, ub_list=ub_list)
433 29085 : DEALLOCATE (map_ub_kind)
434 : END IF
435 : END DO
436 2643 : CALL timestop(handle2)
437 :
438 2643 : END SUBROUTINE force_field_unique_ub
439 :
440 : ! **************************************************************************************************
441 : !> \brief Determine the number of unique torsion kind and allocate torsion_kind_set
442 : !> \param particle_set ...
443 : !> \param molecule_kind_set ...
444 : !> \param molecule_set ...
445 : !> \param ff_type ...
446 : !> \param iw ...
447 : ! **************************************************************************************************
448 2643 : SUBROUTINE force_field_unique_tors(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
449 :
450 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
451 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
452 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
453 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
454 : INTEGER, INTENT(IN) :: iw
455 :
456 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_tors'
457 :
458 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, &
459 : name_atm_b2, name_atm_c, name_atm_c2, &
460 : name_atm_d, name_atm_d2
461 : INTEGER :: atm_a, atm_b, atm_c, atm_d, counter, &
462 : first, handle2, i, j, k, last, natom, &
463 : ntorsion
464 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
465 2643 : INTEGER, POINTER :: map_torsion_kind(:)
466 : LOGICAL :: chk_reverse, found
467 : TYPE(atomic_kind_type), POINTER :: atomic_kind
468 : TYPE(molecule_kind_type), POINTER :: molecule_kind
469 : TYPE(molecule_type), POINTER :: molecule
470 2643 : TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set
471 2643 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
472 :
473 2643 : CALL timeset(routineN, handle2)
474 :
475 2643 : IF (iw > 0) THEN
476 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
477 239 : "FORCEFIELD| Checking for unique torsion terms"
478 : END IF
479 :
480 : ! Now decide whether we need to check D-C-B-A type combination in addtion to usual A-B-C-D
481 : ! We don't need it for Amber FF
482 2643 : chk_reverse = (ff_type%ff_type /= do_ff_amber)
483 :
484 75465 : DO i = 1, SIZE(molecule_kind_set)
485 72822 : molecule_kind => molecule_kind_set(i)
486 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
487 : molecule_list=molecule_list, &
488 : natom=natom, &
489 72822 : ntorsion=ntorsion, torsion_list=torsion_list)
490 72822 : molecule => molecule_set(molecule_list(1))
491 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
492 148287 : IF (ntorsion > 0) THEN
493 16596 : ALLOCATE (map_torsion_kind(ntorsion))
494 5532 : counter = 0
495 5532 : IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN
496 320 : DO j = 1, ntorsion
497 320 : map_torsion_kind(j) = j
498 : END DO
499 8 : counter = ntorsion
500 : ELSE
501 160581 : DO j = 1, ntorsion
502 155057 : atm_a = torsion_list(j)%a
503 155057 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
504 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
505 155057 : name=name_atm_a)
506 155057 : atm_b = torsion_list(j)%b
507 155057 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
508 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
509 155057 : name=name_atm_b)
510 155057 : atm_c = torsion_list(j)%c
511 155057 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
512 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
513 155057 : name=name_atm_c)
514 155057 : atm_d = torsion_list(j)%d
515 155057 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
516 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
517 155057 : name=name_atm_d)
518 155057 : found = .FALSE.
519 2930642 : DO k = 1, j - 1
520 2838283 : atm_a = torsion_list(k)%a
521 2838283 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
522 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
523 2838283 : name=name_atm_a2)
524 2838283 : atm_b = torsion_list(k)%b
525 2838283 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
526 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
527 2838283 : name=name_atm_b2)
528 2838283 : atm_c = torsion_list(k)%c
529 2838283 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
530 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
531 2838283 : name=name_atm_c2)
532 2838283 : atm_d = torsion_list(k)%d
533 2838283 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
534 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
535 2838283 : name=name_atm_d2)
536 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
537 : ((name_atm_b) == (name_atm_b2)) .AND. &
538 : ((name_atm_c) == (name_atm_c2)) .AND. &
539 2838283 : ((name_atm_d) == (name_atm_d2))) .OR. &
540 : (chk_reverse .AND. &
541 : ((name_atm_a) == (name_atm_d2)) .AND. &
542 : ((name_atm_b) == (name_atm_c2)) .AND. &
543 : ((name_atm_c) == (name_atm_b2)) .AND. &
544 92359 : ((name_atm_d) == (name_atm_a2)))) THEN
545 62698 : found = .TRUE.
546 62698 : map_torsion_kind(j) = map_torsion_kind(k)
547 : EXIT
548 : END IF
549 : END DO
550 5524 : IF (.NOT. found) THEN
551 92359 : counter = counter + 1
552 92359 : map_torsion_kind(j) = counter
553 : END IF
554 : END DO
555 : END IF
556 5532 : NULLIFY (torsion_kind_set)
557 5532 : CALL allocate_torsion_kind_set(torsion_kind_set, counter)
558 160901 : DO j = 1, ntorsion
559 160901 : torsion_list(j)%torsion_kind => torsion_kind_set(map_torsion_kind(j))
560 : END DO
561 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
562 5532 : torsion_kind_set=torsion_kind_set, torsion_list=torsion_list)
563 5532 : DEALLOCATE (map_torsion_kind)
564 : END IF
565 : END DO
566 :
567 2643 : CALL timestop(handle2)
568 :
569 2643 : END SUBROUTINE force_field_unique_tors
570 :
571 : ! **************************************************************************************************
572 : !> \brief Determine the number of unique impr kind and allocate impr_kind_set
573 : !> \param particle_set ...
574 : !> \param molecule_kind_set ...
575 : !> \param molecule_set ...
576 : !> \param ff_type ...
577 : !> \param iw ...
578 : ! **************************************************************************************************
579 2643 : SUBROUTINE force_field_unique_impr(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
580 :
581 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
582 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
583 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
584 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
585 : INTEGER, INTENT(IN) :: iw
586 :
587 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_impr'
588 :
589 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, &
590 : name_atm_b2, name_atm_c, name_atm_c2, &
591 : name_atm_d, name_atm_d2
592 : INTEGER :: atm_a, atm_b, atm_c, atm_d, counter, &
593 : first, handle2, i, j, k, last, natom, &
594 : nimpr
595 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
596 2643 : INTEGER, POINTER :: map_impr_kind(:)
597 : LOGICAL :: found
598 : TYPE(atomic_kind_type), POINTER :: atomic_kind
599 2643 : TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set
600 2643 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list
601 : TYPE(molecule_kind_type), POINTER :: molecule_kind
602 : TYPE(molecule_type), POINTER :: molecule
603 :
604 2643 : CALL timeset(routineN, handle2)
605 :
606 2643 : IF (iw > 0) THEN
607 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
608 239 : "FORCEFIELD| Checking for unique improper terms"
609 : END IF
610 :
611 75465 : DO i = 1, SIZE(molecule_kind_set)
612 72822 : molecule_kind => molecule_kind_set(i)
613 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
614 : molecule_list=molecule_list, &
615 : natom=natom, &
616 72822 : nimpr=nimpr, impr_list=impr_list)
617 72822 : molecule => molecule_set(molecule_list(1))
618 :
619 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
620 :
621 148287 : IF (nimpr > 0) THEN
622 5016 : ALLOCATE (map_impr_kind(nimpr))
623 1672 : counter = 0
624 1672 : IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN
625 0 : DO j = 1, nimpr
626 0 : map_impr_kind(j) = j
627 : END DO
628 0 : counter = nimpr
629 : ELSE
630 6984 : DO j = 1, nimpr
631 5312 : atm_a = impr_list(j)%a
632 5312 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
633 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
634 5312 : name=name_atm_a)
635 5312 : atm_b = impr_list(j)%b
636 5312 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
637 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
638 5312 : name=name_atm_b)
639 5312 : atm_c = impr_list(j)%c
640 5312 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
641 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
642 5312 : name=name_atm_c)
643 5312 : atm_d = impr_list(j)%d
644 5312 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
645 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
646 5312 : name=name_atm_d)
647 5312 : found = .FALSE.
648 18542 : DO k = 1, j - 1
649 13834 : atm_a = impr_list(k)%a
650 13834 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
651 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
652 13834 : name=name_atm_a2)
653 13834 : atm_b = impr_list(k)%b
654 13834 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
655 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
656 13834 : name=name_atm_b2)
657 13834 : atm_c = impr_list(k)%c
658 13834 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
659 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
660 13834 : name=name_atm_c2)
661 13834 : atm_d = impr_list(k)%d
662 13834 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
663 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
664 13834 : name=name_atm_d2)
665 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
666 : ((name_atm_b) == (name_atm_b2)) .AND. &
667 : ((name_atm_c) == (name_atm_c2)) .AND. &
668 13834 : ((name_atm_d) == (name_atm_d2))) .OR. &
669 : (((name_atm_a) == (name_atm_a2)) .AND. &
670 : ((name_atm_b) == (name_atm_b2)) .AND. &
671 : ((name_atm_c) == (name_atm_d2)) .AND. &
672 4708 : ((name_atm_d) == (name_atm_c2)))) THEN
673 604 : found = .TRUE.
674 604 : map_impr_kind(j) = map_impr_kind(k)
675 : EXIT
676 : END IF
677 : END DO
678 1672 : IF (.NOT. found) THEN
679 4708 : counter = counter + 1
680 4708 : map_impr_kind(j) = counter
681 : END IF
682 : END DO
683 : END IF
684 1672 : NULLIFY (impr_kind_set)
685 1672 : CALL allocate_impr_kind_set(impr_kind_set, counter)
686 6984 : DO j = 1, nimpr
687 6984 : impr_list(j)%impr_kind => impr_kind_set(map_impr_kind(j))
688 : END DO
689 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
690 1672 : impr_kind_set=impr_kind_set, impr_list=impr_list)
691 1672 : DEALLOCATE (map_impr_kind)
692 : END IF
693 : END DO
694 2643 : CALL timestop(handle2)
695 :
696 2643 : END SUBROUTINE force_field_unique_impr
697 :
698 : ! **************************************************************************************************
699 : !> \brief Determine the number of unique opbend kind and allocate opbend_kind_set
700 : !> based on the present impropers. With each improper, there also
701 : !> corresponds a opbend
702 : !> \param particle_set ...
703 : !> \param molecule_kind_set ...
704 : !> \param molecule_set ...
705 : !> \param ff_type ...
706 : !> \param iw ...
707 : ! **************************************************************************************************
708 2643 : SUBROUTINE force_field_unique_opbend(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
709 :
710 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
711 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
712 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
713 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
714 : INTEGER, INTENT(IN) :: iw
715 :
716 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_opbend'
717 :
718 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, &
719 : name_atm_b2, name_atm_c, name_atm_c2, &
720 : name_atm_d, name_atm_d2
721 : INTEGER :: atm_a, atm_b, atm_c, atm_d, counter, &
722 : first, handle2, i, j, k, last, natom, &
723 : nopbend
724 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
725 2643 : INTEGER, POINTER :: map_opbend_kind(:)
726 : LOGICAL :: found
727 : TYPE(atomic_kind_type), POINTER :: atomic_kind
728 : TYPE(molecule_kind_type), POINTER :: molecule_kind
729 : TYPE(molecule_type), POINTER :: molecule
730 2643 : TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set
731 2643 : TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list
732 :
733 2643 : CALL timeset(routineN, handle2)
734 :
735 2643 : IF (iw > 0) THEN
736 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
737 239 : "FORCEFIELD| Checking for unique out-of-plane bend terms"
738 : END IF
739 :
740 75465 : DO i = 1, SIZE(molecule_kind_set)
741 72822 : molecule_kind => molecule_kind_set(i)
742 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
743 : molecule_list=molecule_list, &
744 : natom=natom, &
745 72822 : nopbend=nopbend, opbend_list=opbend_list)
746 72822 : molecule => molecule_set(molecule_list(1))
747 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
748 148287 : IF (nopbend > 0) THEN
749 5016 : ALLOCATE (map_opbend_kind(nopbend))
750 1672 : counter = 0
751 1672 : IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN
752 0 : DO j = 1, nopbend
753 0 : map_opbend_kind(j) = j
754 : END DO
755 0 : counter = nopbend
756 : ELSE
757 6984 : DO j = 1, nopbend
758 5312 : atm_a = opbend_list(j)%a
759 5312 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
760 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
761 5312 : name=name_atm_a)
762 5312 : atm_b = opbend_list(j)%b
763 5312 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
764 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
765 5312 : name=name_atm_b)
766 5312 : atm_c = opbend_list(j)%c
767 5312 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
768 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
769 5312 : name=name_atm_c)
770 5312 : atm_d = opbend_list(j)%d
771 5312 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
772 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
773 5312 : name=name_atm_d)
774 5312 : found = .FALSE.
775 18542 : DO k = 1, j - 1
776 13834 : atm_a = opbend_list(k)%a
777 13834 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
778 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
779 13834 : name=name_atm_a2)
780 13834 : atm_b = opbend_list(k)%b
781 13834 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
782 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
783 13834 : name=name_atm_b2)
784 13834 : atm_c = opbend_list(k)%c
785 13834 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
786 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
787 13834 : name=name_atm_c2)
788 13834 : atm_d = opbend_list(k)%d
789 13834 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
790 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
791 13834 : name=name_atm_d2)
792 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
793 : ((name_atm_b) == (name_atm_b2)) .AND. &
794 : ((name_atm_c) == (name_atm_c2)) .AND. &
795 13834 : ((name_atm_d) == (name_atm_d2))) .OR. &
796 : (((name_atm_a) == (name_atm_a2)) .AND. &
797 : ((name_atm_b) == (name_atm_c2)) .AND. &
798 : ((name_atm_c) == (name_atm_b2)) .AND. &
799 4708 : ((name_atm_d) == (name_atm_d2)))) THEN
800 604 : found = .TRUE.
801 604 : map_opbend_kind(j) = map_opbend_kind(k)
802 : EXIT
803 : END IF
804 : END DO
805 1672 : IF (.NOT. found) THEN
806 4708 : counter = counter + 1
807 4708 : map_opbend_kind(j) = counter
808 : END IF
809 : END DO
810 : END IF
811 1672 : NULLIFY (opbend_kind_set)
812 1672 : CALL allocate_opbend_kind_set(opbend_kind_set, counter)
813 6984 : DO j = 1, nopbend
814 6984 : opbend_list(j)%opbend_kind => opbend_kind_set(map_opbend_kind(j))
815 : END DO
816 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
817 1672 : opbend_kind_set=opbend_kind_set, opbend_list=opbend_list)
818 1672 : DEALLOCATE (map_opbend_kind)
819 : END IF
820 : END DO
821 2643 : CALL timestop(handle2)
822 :
823 2643 : END SUBROUTINE force_field_unique_opbend
824 :
825 : ! **************************************************************************************************
826 : !> \brief Pack in bonds information needed for the force_field
827 : !> \param particle_set ...
828 : !> \param molecule_kind_set ...
829 : !> \param molecule_set ...
830 : !> \param fatal ...
831 : !> \param Ainfo ...
832 : !> \param chm_info ...
833 : !> \param inp_info ...
834 : !> \param gro_info ...
835 : !> \param amb_info ...
836 : !> \param iw ...
837 : ! **************************************************************************************************
838 2643 : SUBROUTINE force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, fatal, Ainfo, &
839 : chm_info, inp_info, gro_info, amb_info, iw)
840 :
841 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
842 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
843 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
844 : LOGICAL :: fatal
845 : CHARACTER(LEN=default_string_length), &
846 : DIMENSION(:), POINTER :: Ainfo
847 : TYPE(charmm_info_type), POINTER :: chm_info
848 : TYPE(input_info_type), POINTER :: inp_info
849 : TYPE(gromos_info_type), POINTER :: gro_info
850 : TYPE(amber_info_type), POINTER :: amb_info
851 : INTEGER, INTENT(IN) :: iw
852 :
853 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bond'
854 :
855 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b
856 : INTEGER :: atm_a, atm_b, first, handle2, i, itype, &
857 : j, k, last, natom, nbond
858 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
859 : LOGICAL :: found, only_qm
860 : TYPE(atomic_kind_type), POINTER :: atomic_kind
861 2643 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
862 : TYPE(molecule_kind_type), POINTER :: molecule_kind
863 : TYPE(molecule_type), POINTER :: molecule
864 :
865 2643 : CALL timeset(routineN, handle2)
866 :
867 2643 : IF (iw > 0) THEN
868 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
869 239 : "FORCEFIELD| Checking for bond terms"
870 : END IF
871 :
872 75465 : DO i = 1, SIZE(molecule_kind_set)
873 72822 : molecule_kind => molecule_kind_set(i)
874 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
875 : molecule_list=molecule_list, &
876 : natom=natom, &
877 72822 : nbond=nbond, bond_list=bond_list)
878 72822 : molecule => molecule_set(molecule_list(1))
879 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
880 186997 : DO j = 1, nbond
881 114175 : atm_a = bond_list(j)%a
882 114175 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
883 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
884 114175 : name=name_atm_a)
885 114175 : atm_b = bond_list(j)%b
886 114175 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
887 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
888 114175 : name=name_atm_b)
889 114175 : found = .FALSE.
890 114175 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
891 114175 : CALL uppercase(name_atm_a)
892 114175 : CALL uppercase(name_atm_b)
893 :
894 : ! loop over params from GROMOS
895 114175 : IF (ASSOCIATED(gro_info%bond_k)) THEN
896 128 : k = SIZE(gro_info%bond_k)
897 128 : itype = bond_list(j)%itype
898 128 : IF (itype <= k) THEN
899 104 : bond_list(j)%bond_kind%k(1) = gro_info%bond_k(itype)
900 104 : bond_list(j)%bond_kind%r0 = gro_info%bond_r0(itype)
901 : ELSE
902 24 : itype = itype - k
903 24 : bond_list(j)%bond_kind%k(1) = gro_info%solvent_k(itype)
904 24 : bond_list(j)%bond_kind%r0 = gro_info%solvent_r0(itype)
905 : END IF
906 128 : bond_list(j)%bond_kind%id_type = gro_info%ff_gromos_type
907 128 : bond_list(j)%id_type = gro_info%ff_gromos_type
908 128 : found = .TRUE.
909 : END IF
910 :
911 : ! loop over params from CHARMM
912 114175 : IF (ASSOCIATED(chm_info%bond_a)) THEN
913 1449348 : DO k = 1, SIZE(chm_info%bond_a)
914 : IF ((((chm_info%bond_a(k)) == (name_atm_a)) .AND. &
915 1449324 : ((chm_info%bond_b(k)) == (name_atm_b))) .OR. &
916 : (((chm_info%bond_a(k)) == (name_atm_b)) .AND. &
917 24 : ((chm_info%bond_b(k)) == (name_atm_a)))) THEN
918 41447 : bond_list(j)%bond_kind%id_type = do_ff_charmm
919 41447 : bond_list(j)%bond_kind%k(1) = chm_info%bond_k(k)
920 41447 : bond_list(j)%bond_kind%r0 = chm_info%bond_r0(k)
921 41447 : CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b)
922 41447 : found = .TRUE.
923 41447 : EXIT
924 : END IF
925 : END DO
926 : END IF
927 :
928 : ! loop over params from AMBER
929 114175 : IF (ASSOCIATED(amb_info%bond_a)) THEN
930 5716862 : DO k = 1, SIZE(amb_info%bond_a)
931 : IF ((((amb_info%bond_a(k)) == (name_atm_a)) .AND. &
932 5716862 : ((amb_info%bond_b(k)) == (name_atm_b))) .OR. &
933 : (((amb_info%bond_a(k)) == (name_atm_b)) .AND. &
934 0 : ((amb_info%bond_b(k)) == (name_atm_a)))) THEN
935 64808 : bond_list(j)%bond_kind%id_type = do_ff_amber
936 64808 : bond_list(j)%bond_kind%k(1) = amb_info%bond_k(k)
937 64808 : bond_list(j)%bond_kind%r0 = amb_info%bond_r0(k)
938 64808 : CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b)
939 64808 : found = .TRUE.
940 64808 : EXIT
941 : END IF
942 : END DO
943 : END IF
944 :
945 : ! always have the input param last to overwrite all the other ones
946 114175 : IF (ASSOCIATED(inp_info%bond_a)) THEN
947 9672 : DO k = 1, SIZE(inp_info%bond_a)
948 : IF ((((inp_info%bond_a(k)) == (name_atm_a)) .AND. &
949 9626 : ((inp_info%bond_b(k)) == (name_atm_b))) .OR. &
950 : (((inp_info%bond_a(k)) == (name_atm_b)) .AND. &
951 46 : ((inp_info%bond_b(k)) == (name_atm_a)))) THEN
952 7800 : bond_list(j)%bond_kind%id_type = inp_info%bond_kind(k)
953 54600 : bond_list(j)%bond_kind%k(:) = inp_info%bond_k(:, k)
954 7800 : bond_list(j)%bond_kind%r0 = inp_info%bond_r0(k)
955 7800 : bond_list(j)%bond_kind%cs = inp_info%bond_cs(k)
956 7800 : CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b)
957 7800 : found = .TRUE.
958 7800 : EXIT
959 : END IF
960 : END DO
961 : END IF
962 :
963 114175 : IF (.NOT. found) CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
964 : atm2=TRIM(name_atm_b), &
965 : fatal=fatal, &
966 : type_name="Bond", &
967 16 : array=Ainfo)
968 : ! QM/MM modifications
969 186997 : IF (only_qm) THEN
970 2082 : bond_list(j)%id_type = do_ff_undef
971 2082 : bond_list(j)%bond_kind%id_type = do_ff_undef
972 : END IF
973 : END DO
974 :
975 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
976 148287 : bond_list=bond_list)
977 :
978 : END DO
979 :
980 2643 : CALL timestop(handle2)
981 :
982 2643 : END SUBROUTINE force_field_pack_bond
983 :
984 : ! **************************************************************************************************
985 : !> \brief Pack in bends information needed for the force_field
986 : !> \param particle_set ...
987 : !> \param molecule_kind_set ...
988 : !> \param molecule_set ...
989 : !> \param fatal ...
990 : !> \param Ainfo ...
991 : !> \param chm_info ...
992 : !> \param inp_info ...
993 : !> \param gro_info ...
994 : !> \param amb_info ...
995 : !> \param iw ...
996 : ! **************************************************************************************************
997 2643 : SUBROUTINE force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, fatal, Ainfo, &
998 : chm_info, inp_info, gro_info, amb_info, iw)
999 :
1000 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1001 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1002 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1003 : LOGICAL :: fatal
1004 : CHARACTER(LEN=default_string_length), &
1005 : DIMENSION(:), POINTER :: Ainfo
1006 : TYPE(charmm_info_type), POINTER :: chm_info
1007 : TYPE(input_info_type), POINTER :: inp_info
1008 : TYPE(gromos_info_type), POINTER :: gro_info
1009 : TYPE(amber_info_type), POINTER :: amb_info
1010 : INTEGER, INTENT(IN) :: iw
1011 :
1012 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bend'
1013 :
1014 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c
1015 : INTEGER :: atm_a, atm_b, atm_c, first, handle2, i, &
1016 : itype, j, k, l, last, natom, nbend
1017 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1018 : LOGICAL :: found, only_qm
1019 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1020 2643 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
1021 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1022 : TYPE(molecule_type), POINTER :: molecule
1023 :
1024 2643 : CALL timeset(routineN, handle2)
1025 :
1026 2643 : IF (iw > 0) THEN
1027 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1028 239 : "FORCEFIELD| Checking for bend terms"
1029 : END IF
1030 :
1031 75465 : DO i = 1, SIZE(molecule_kind_set)
1032 72822 : molecule_kind => molecule_kind_set(i)
1033 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1034 : molecule_list=molecule_list, &
1035 : natom=natom, &
1036 72822 : nbend=nbend, bend_list=bend_list)
1037 72822 : molecule => molecule_set(molecule_list(1))
1038 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1039 211916 : DO j = 1, nbend
1040 139094 : atm_a = bend_list(j)%a
1041 139094 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1042 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1043 139094 : name=name_atm_a)
1044 139094 : atm_b = bend_list(j)%b
1045 139094 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1046 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1047 139094 : name=name_atm_b)
1048 139094 : atm_c = bend_list(j)%c
1049 139094 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1050 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1051 139094 : name=name_atm_c)
1052 139094 : found = .FALSE.
1053 139094 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c)
1054 139094 : CALL uppercase(name_atm_a)
1055 139094 : CALL uppercase(name_atm_b)
1056 139094 : CALL uppercase(name_atm_c)
1057 :
1058 : ! loop over params from GROMOS
1059 139094 : IF (ASSOCIATED(gro_info%bend_k)) THEN
1060 156 : k = SIZE(gro_info%bend_k)
1061 156 : itype = bend_list(j)%itype
1062 156 : IF (itype > 0) THEN
1063 156 : bend_list(j)%bend_kind%k = gro_info%bend_k(itype)
1064 156 : bend_list(j)%bend_kind%theta0 = gro_info%bend_theta0(itype)
1065 : ELSE
1066 0 : bend_list(j)%bend_kind%k = gro_info%bend_k(itype/k)
1067 0 : bend_list(j)%bend_kind%theta0 = gro_info%bend_theta0(itype/k)
1068 : END IF
1069 156 : bend_list(j)%bend_kind%id_type = gro_info%ff_gromos_type
1070 156 : bend_list(j)%id_type = gro_info%ff_gromos_type
1071 156 : found = .TRUE.
1072 : END IF
1073 :
1074 : ! loop over params from CHARMM
1075 139094 : IF (ASSOCIATED(chm_info%bend_a)) THEN
1076 6045171 : DO k = 1, SIZE(chm_info%bend_a)
1077 : IF ((((chm_info%bend_a(k)) == (name_atm_a)) .AND. &
1078 : ((chm_info%bend_b(k)) == (name_atm_b)) .AND. &
1079 6045097 : ((chm_info%bend_c(k)) == (name_atm_c))) .OR. &
1080 : (((chm_info%bend_a(k)) == (name_atm_c)) .AND. &
1081 : ((chm_info%bend_b(k)) == (name_atm_b)) .AND. &
1082 74 : ((chm_info%bend_c(k)) == (name_atm_a)))) THEN
1083 67523 : bend_list(j)%bend_kind%id_type = do_ff_charmm
1084 67523 : bend_list(j)%bend_kind%k = chm_info%bend_k(k)
1085 67523 : bend_list(j)%bend_kind%theta0 = chm_info%bend_theta0(k)
1086 : CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, &
1087 67523 : name_atm_c)
1088 67523 : found = .TRUE.
1089 67523 : EXIT
1090 : END IF
1091 : END DO
1092 : END IF
1093 :
1094 : ! loop over params from AMBER
1095 139094 : IF (ASSOCIATED(amb_info%bend_a)) THEN
1096 10981138 : DO k = 1, SIZE(amb_info%bend_a)
1097 : IF ((((amb_info%bend_a(k)) == (name_atm_a)) .AND. &
1098 : ((amb_info%bend_b(k)) == (name_atm_b)) .AND. &
1099 10981138 : ((amb_info%bend_c(k)) == (name_atm_c))) .OR. &
1100 : (((amb_info%bend_a(k)) == (name_atm_c)) .AND. &
1101 : ((amb_info%bend_b(k)) == (name_atm_b)) .AND. &
1102 0 : ((amb_info%bend_c(k)) == (name_atm_a)))) THEN
1103 59540 : bend_list(j)%bend_kind%id_type = do_ff_amber
1104 59540 : bend_list(j)%bend_kind%k = amb_info%bend_k(k)
1105 59540 : bend_list(j)%bend_kind%theta0 = amb_info%bend_theta0(k)
1106 : CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, &
1107 59540 : name_atm_c)
1108 59540 : found = .TRUE.
1109 59540 : EXIT
1110 : END IF
1111 : END DO
1112 : END IF
1113 :
1114 : ! always have the input param last to overwrite all the other ones
1115 139094 : IF (ASSOCIATED(inp_info%bend_a)) THEN
1116 25741 : DO k = 1, SIZE(inp_info%bend_a)
1117 : IF ((((inp_info%bend_a(k)) == (name_atm_a)) .AND. &
1118 : ((inp_info%bend_b(k)) == (name_atm_b)) .AND. &
1119 25725 : ((inp_info%bend_c(k)) == (name_atm_c))) .OR. &
1120 : (((inp_info%bend_a(k)) == (name_atm_c)) .AND. &
1121 : ((inp_info%bend_b(k)) == (name_atm_b)) .AND. &
1122 16 : ((inp_info%bend_c(k)) == (name_atm_a)))) THEN
1123 11875 : bend_list(j)%bend_kind%id_type = inp_info%bend_kind(k)
1124 11875 : bend_list(j)%bend_kind%k = inp_info%bend_k(k)
1125 11875 : bend_list(j)%bend_kind%theta0 = inp_info%bend_theta0(k)
1126 11875 : bend_list(j)%bend_kind%cb = inp_info%bend_cb(k)
1127 11875 : bend_list(j)%bend_kind%r012 = inp_info%bend_r012(k)
1128 11875 : bend_list(j)%bend_kind%r032 = inp_info%bend_r032(k)
1129 11875 : bend_list(j)%bend_kind%kbs12 = inp_info%bend_kbs12(k)
1130 11875 : bend_list(j)%bend_kind%kbs32 = inp_info%bend_kbs32(k)
1131 11875 : bend_list(j)%bend_kind%kss = inp_info%bend_kss(k)
1132 11875 : bend_list(j)%bend_kind%legendre%order = inp_info%bend_legendre(k)%order
1133 11875 : IF (bend_list(j)%bend_kind%legendre%order /= 0) THEN
1134 11875 : IF (ASSOCIATED(bend_list(j)%bend_kind%legendre%coeffs)) THEN
1135 9554 : DEALLOCATE (bend_list(j)%bend_kind%legendre%coeffs)
1136 : END IF
1137 35625 : ALLOCATE (bend_list(j)%bend_kind%legendre%coeffs(bend_list(j)%bend_kind%legendre%order))
1138 23990 : DO l = 1, bend_list(j)%bend_kind%legendre%order
1139 23990 : bend_list(j)%bend_kind%legendre%coeffs(l) = inp_info%bend_legendre(k)%coeffs(l)
1140 : END DO
1141 : END IF
1142 : CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, &
1143 11875 : name_atm_c)
1144 11875 : found = .TRUE.
1145 11875 : EXIT
1146 : END IF
1147 : END DO
1148 : END IF
1149 :
1150 139094 : IF (.NOT. found) CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
1151 : atm2=TRIM(name_atm_b), &
1152 : atm3=TRIM(name_atm_c), &
1153 : fatal=fatal, &
1154 : type_name="Angle", &
1155 8 : array=Ainfo)
1156 : ! QM/MM modifications
1157 211916 : IF (only_qm) THEN
1158 1918 : bend_list(j)%id_type = do_ff_undef
1159 1918 : bend_list(j)%bend_kind%id_type = do_ff_undef
1160 : END IF
1161 : END DO
1162 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1163 148287 : bend_list=bend_list)
1164 : END DO
1165 2643 : CALL timestop(handle2)
1166 :
1167 2643 : END SUBROUTINE force_field_pack_bend
1168 :
1169 : ! **************************************************************************************************
1170 : !> \brief Pack in Urey-Bradley information needed for the force_field
1171 : !> \param particle_set ...
1172 : !> \param molecule_kind_set ...
1173 : !> \param molecule_set ...
1174 : !> \param Ainfo ...
1175 : !> \param chm_info ...
1176 : !> \param inp_info ...
1177 : !> \param iw ...
1178 : ! **************************************************************************************************
1179 2643 : SUBROUTINE force_field_pack_ub(particle_set, molecule_kind_set, molecule_set, &
1180 : Ainfo, chm_info, inp_info, iw)
1181 :
1182 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1183 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1184 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1185 : CHARACTER(LEN=default_string_length), &
1186 : DIMENSION(:), POINTER :: Ainfo
1187 : TYPE(charmm_info_type), POINTER :: chm_info
1188 : TYPE(input_info_type), POINTER :: inp_info
1189 : INTEGER :: iw
1190 :
1191 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_ub'
1192 :
1193 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c
1194 : INTEGER :: atm_a, atm_b, atm_c, first, handle2, i, &
1195 : j, k, last, natom, nub
1196 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1197 : LOGICAL :: found, only_qm
1198 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1199 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1200 : TYPE(molecule_type), POINTER :: molecule
1201 2643 : TYPE(ub_type), DIMENSION(:), POINTER :: ub_list
1202 :
1203 2643 : CALL timeset(routineN, handle2)
1204 :
1205 2643 : IF (iw > 0) THEN
1206 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1207 239 : "FORCEFIELD| Checking for Urey-Bradley (UB) terms"
1208 : END IF
1209 :
1210 75465 : DO i = 1, SIZE(molecule_kind_set)
1211 72822 : molecule_kind => molecule_kind_set(i)
1212 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1213 : molecule_list=molecule_list, &
1214 : natom=natom, &
1215 72822 : nub=nub, ub_list=ub_list)
1216 72822 : molecule => molecule_set(molecule_list(1))
1217 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1218 211758 : DO j = 1, nub
1219 138936 : atm_a = ub_list(j)%a
1220 138936 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1221 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1222 138936 : name=name_atm_a)
1223 138936 : atm_b = ub_list(j)%b
1224 138936 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1225 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1226 138936 : name=name_atm_b)
1227 138936 : atm_c = ub_list(j)%c
1228 138936 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1229 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1230 138936 : name=name_atm_c)
1231 138936 : found = .FALSE.
1232 138936 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c)
1233 138936 : CALL uppercase(name_atm_a)
1234 138936 : CALL uppercase(name_atm_b)
1235 138936 : CALL uppercase(name_atm_c)
1236 :
1237 : ! Loop over params from GROMOS
1238 : ! ikuo - None that I know...
1239 :
1240 : ! Loop over params from CHARMM
1241 138936 : IF (ASSOCIATED(chm_info%ub_a)) THEN
1242 3842528 : DO k = 1, SIZE(chm_info%ub_a)
1243 : IF ((((chm_info%ub_a(k)) == (name_atm_a)) .AND. &
1244 : ((chm_info%ub_b(k)) == (name_atm_b)) .AND. &
1245 3818446 : ((chm_info%ub_c(k)) == (name_atm_c))) .OR. &
1246 : (((chm_info%ub_a(k)) == (name_atm_c)) .AND. &
1247 : ((chm_info%ub_b(k)) == (name_atm_b)) .AND. &
1248 24082 : ((chm_info%ub_c(k)) == (name_atm_a)))) THEN
1249 20692 : ub_list(j)%ub_kind%id_type = do_ff_charmm
1250 20692 : ub_list(j)%ub_kind%k(1) = chm_info%ub_k(k)
1251 20692 : ub_list(j)%ub_kind%r0 = chm_info%ub_r0(k)
1252 20692 : IF (iw > 0) THEN
1253 : WRITE (UNIT=iw, FMT="(T2,A)") &
1254 : "FORCEFIELD| Found Urey-Bradley term (CHARMM) for the atomic kinds "// &
1255 138 : TRIM(name_atm_a)//", "//TRIM(name_atm_b)//" and "//TRIM(name_atm_c)
1256 : END IF
1257 : CALL issue_duplications(found, "Urey-Bradley", name_atm_a, &
1258 20692 : name_atm_b, name_atm_c)
1259 20692 : found = .TRUE.
1260 20692 : EXIT
1261 : END IF
1262 : END DO
1263 : END IF
1264 :
1265 : ! Loop over params from AMBER
1266 : ! teo - None that I know...
1267 :
1268 : ! Always have the input param last to overwrite all the other ones
1269 138936 : IF (ASSOCIATED(inp_info%ub_a)) THEN
1270 45592 : DO k = 1, SIZE(inp_info%ub_a)
1271 : IF ((((inp_info%ub_a(k)) == (name_atm_a)) .AND. &
1272 : ((inp_info%ub_b(k)) == (name_atm_b)) .AND. &
1273 33709 : ((inp_info%ub_c(k)) == (name_atm_c))) .OR. &
1274 : (((inp_info%ub_a(k)) == (name_atm_c)) .AND. &
1275 : ((inp_info%ub_b(k)) == (name_atm_b)) .AND. &
1276 11883 : ((inp_info%ub_c(k)) == (name_atm_a)))) THEN
1277 8 : ub_list(j)%ub_kind%id_type = inp_info%ub_kind(k)
1278 56 : ub_list(j)%ub_kind%k(:) = inp_info%ub_k(:, k)
1279 8 : ub_list(j)%ub_kind%r0 = inp_info%ub_r0(k)
1280 8 : IF (iw > 0) THEN
1281 : WRITE (UNIT=iw, FMT="(T2,A)") &
1282 : "FORCEFIELD| Found Urey-Bradley term (input) for the atomic kinds "// &
1283 0 : TRIM(name_atm_a)//", "//TRIM(name_atm_b)//" and "//TRIM(name_atm_c)
1284 : END IF
1285 : CALL issue_duplications(found, "Urey-Bradley", name_atm_a, &
1286 8 : name_atm_b, name_atm_c)
1287 8 : found = .TRUE.
1288 8 : EXIT
1289 : END IF
1290 : END DO
1291 : END IF
1292 :
1293 138936 : IF (.NOT. found) THEN
1294 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
1295 : atm2=TRIM(name_atm_b), &
1296 : atm3=TRIM(name_atm_c), &
1297 : type_name="Urey-Bradley", &
1298 118236 : array=Ainfo)
1299 118236 : ub_list(j)%id_type = do_ff_undef
1300 118236 : ub_list(j)%ub_kind%id_type = do_ff_undef
1301 472944 : ub_list(j)%ub_kind%k = 0.0_dp
1302 118236 : ub_list(j)%ub_kind%r0 = 0.0_dp
1303 : END IF
1304 :
1305 : ! QM/MM modifications
1306 211758 : IF (only_qm) THEN
1307 1918 : ub_list(j)%id_type = do_ff_undef
1308 1918 : ub_list(j)%ub_kind%id_type = do_ff_undef
1309 : END IF
1310 : END DO
1311 :
1312 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1313 148287 : ub_list=ub_list)
1314 :
1315 : END DO
1316 :
1317 2643 : CALL timestop(handle2)
1318 :
1319 2643 : END SUBROUTINE force_field_pack_ub
1320 :
1321 : ! **************************************************************************************************
1322 : !> \brief Pack in torsion information needed for the force_field
1323 : !> \param particle_set ...
1324 : !> \param molecule_kind_set ...
1325 : !> \param molecule_set ...
1326 : !> \param Ainfo ...
1327 : !> \param chm_info ...
1328 : !> \param inp_info ...
1329 : !> \param gro_info ...
1330 : !> \param amb_info ...
1331 : !> \param iw ...
1332 : ! **************************************************************************************************
1333 2643 : SUBROUTINE force_field_pack_tors(particle_set, molecule_kind_set, molecule_set, &
1334 : Ainfo, chm_info, inp_info, gro_info, amb_info, iw)
1335 :
1336 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1337 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1338 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1339 : CHARACTER(LEN=default_string_length), &
1340 : DIMENSION(:), POINTER :: Ainfo
1341 : TYPE(charmm_info_type), POINTER :: chm_info
1342 : TYPE(input_info_type), POINTER :: inp_info
1343 : TYPE(gromos_info_type), POINTER :: gro_info
1344 : TYPE(amber_info_type), POINTER :: amb_info
1345 : INTEGER, INTENT(IN) :: iw
1346 :
1347 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_tors'
1348 :
1349 : CHARACTER(LEN=default_string_length) :: ldum, molecule_name, name_atm_a, &
1350 : name_atm_b, name_atm_c, name_atm_d
1351 : INTEGER :: atm_a, atm_b, atm_c, atm_d, first, &
1352 : handle2, i, imul, itype, j, k, k_end, &
1353 : k_start, last, natom, ntorsion, &
1354 : raw_parm_id
1355 : INTEGER, DIMENSION(4) :: glob_atm_id
1356 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1357 : LOGICAL :: found, only_qm
1358 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1359 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1360 : TYPE(molecule_type), POINTER :: molecule
1361 2643 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
1362 :
1363 2643 : CALL timeset(routineN, handle2)
1364 :
1365 2643 : IF (iw > 0) THEN
1366 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1367 239 : "FORCEFIELD| Checking for torsion terms"
1368 : END IF
1369 :
1370 75465 : DO i = 1, SIZE(molecule_kind_set)
1371 72822 : molecule_kind => molecule_kind_set(i)
1372 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1373 : molecule_list=molecule_list, &
1374 : name=molecule_name, &
1375 : natom=natom, &
1376 : ntorsion=ntorsion, &
1377 72822 : torsion_list=torsion_list)
1378 72822 : molecule => molecule_set(molecule_list(1))
1379 : CALL get_molecule(molecule=molecule, &
1380 : first_atom=first, &
1381 72822 : last_atom=last)
1382 228191 : DO j = 1, ntorsion
1383 228191 : IF (torsion_list(j)%torsion_kind%id_type == do_ff_undef) THEN
1384 113867 : atm_a = torsion_list(j)%a
1385 113867 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1386 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1387 113867 : name=name_atm_a)
1388 113867 : atm_b = torsion_list(j)%b
1389 113867 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1390 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1391 113867 : name=name_atm_b)
1392 113867 : atm_c = torsion_list(j)%c
1393 113867 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1394 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1395 113867 : name=name_atm_c)
1396 113867 : atm_d = torsion_list(j)%d
1397 113867 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
1398 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1399 113867 : name=name_atm_d)
1400 113867 : found = .FALSE.
1401 113867 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d)
1402 113867 : CALL uppercase(name_atm_a)
1403 113867 : CALL uppercase(name_atm_b)
1404 113867 : CALL uppercase(name_atm_c)
1405 113867 : CALL uppercase(name_atm_d)
1406 :
1407 : ! Loop over params from GROMOS
1408 113867 : IF (ASSOCIATED(gro_info%torsion_k)) THEN
1409 312 : k = SIZE(gro_info%torsion_k)
1410 312 : itype = torsion_list(j)%itype
1411 312 : IF (itype > 0) THEN
1412 312 : CALL reallocate(torsion_list(j)%torsion_kind%k, 1, 1)
1413 312 : CALL reallocate(torsion_list(j)%torsion_kind%m, 1, 1)
1414 312 : CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, 1)
1415 312 : torsion_list(j)%torsion_kind%nmul = 1
1416 312 : torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype)
1417 312 : torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype)
1418 312 : torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype)
1419 : ELSE
1420 0 : CALL reallocate(torsion_list(j)%torsion_kind%k, 1, 1)
1421 0 : CALL reallocate(torsion_list(j)%torsion_kind%m, 1, 1)
1422 0 : CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, 1)
1423 0 : torsion_list(j)%torsion_kind%nmul = 1
1424 0 : torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype/k)
1425 0 : torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype/k)
1426 0 : torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype/k)
1427 : END IF
1428 312 : torsion_list(j)%torsion_kind%id_type = gro_info%ff_gromos_type
1429 312 : torsion_list(j)%id_type = gro_info%ff_gromos_type
1430 312 : found = .TRUE.
1431 312 : imul = torsion_list(j)%torsion_kind%nmul
1432 : END IF
1433 :
1434 : ! Loop over params from CHARMM
1435 113867 : IF (ASSOCIATED(chm_info%torsion_a)) THEN
1436 20328202 : DO k = 1, SIZE(chm_info%torsion_a)
1437 : IF ((((chm_info%torsion_a(k)) == (name_atm_a)) .AND. &
1438 : ((chm_info%torsion_b(k)) == (name_atm_b)) .AND. &
1439 : ((chm_info%torsion_c(k)) == (name_atm_c)) .AND. &
1440 20273793 : ((chm_info%torsion_d(k)) == (name_atm_d))) .OR. &
1441 : (((chm_info%torsion_a(k)) == (name_atm_d)) .AND. &
1442 : ((chm_info%torsion_b(k)) == (name_atm_c)) .AND. &
1443 : ((chm_info%torsion_c(k)) == (name_atm_b)) .AND. &
1444 54409 : ((chm_info%torsion_d(k)) == (name_atm_a)))) THEN
1445 44224 : imul = torsion_list(j)%torsion_kind%nmul + 1
1446 44224 : CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul)
1447 44224 : CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul)
1448 44224 : CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul)
1449 44224 : torsion_list(j)%torsion_kind%id_type = do_ff_charmm
1450 44224 : torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k)
1451 44224 : torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k)
1452 44224 : torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k)
1453 44224 : torsion_list(j)%torsion_kind%nmul = imul
1454 44224 : found = .TRUE.
1455 : END IF
1456 : END DO
1457 :
1458 54409 : IF (.NOT. found) THEN
1459 6901506 : DO k = 1, SIZE(chm_info%torsion_a)
1460 : IF ((((chm_info%torsion_a(k)) == ("X")) .AND. &
1461 : ((chm_info%torsion_b(k)) == (name_atm_b)) .AND. &
1462 : ((chm_info%torsion_c(k)) == (name_atm_c)) .AND. &
1463 6886624 : ((chm_info%torsion_d(k)) == ("X"))) .OR. &
1464 : (((chm_info%torsion_a(k)) == ("X")) .AND. &
1465 : ((chm_info%torsion_b(k)) == (name_atm_c)) .AND. &
1466 : ((chm_info%torsion_c(k)) == (name_atm_b)) .AND. &
1467 14882 : ((chm_info%torsion_d(k)) == ("X")))) THEN
1468 12990 : imul = torsion_list(j)%torsion_kind%nmul + 1
1469 12990 : CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul)
1470 12990 : CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul)
1471 12990 : CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul)
1472 12990 : torsion_list(j)%torsion_kind%id_type = do_ff_charmm
1473 12990 : torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k)
1474 12990 : torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k)
1475 12990 : torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k)
1476 12990 : torsion_list(j)%torsion_kind%nmul = imul
1477 12990 : found = .TRUE.
1478 : END IF
1479 : END DO
1480 : END IF
1481 : END IF
1482 :
1483 : ! Loop over params from AMBER
1484 : ! Assign real parameters from Amber PRMTOP file using global atom indices
1485 : ! Type-based assignment is prone to errors
1486 113867 : IF (ASSOCIATED(amb_info%torsion_a)) THEN
1487 : ! Get global atom indices
1488 45098 : glob_atm_id(1) = atm_a + first - 1
1489 45098 : glob_atm_id(2) = atm_b + first - 1
1490 45098 : glob_atm_id(3) = atm_c + first - 1
1491 45098 : glob_atm_id(4) = atm_d + first - 1
1492 :
1493 : ! Search sorted array of raw torsion parameters
1494 : ! The array can be too long for linear lookup
1495 : ! Use binary search for first atom index
1496 45098 : k_start = bsearch_leftmost_2d(amb_info%raw_torsion_id, glob_atm_id(1))
1497 45098 : k_end = UBOUND(amb_info%raw_torsion_id, DIM=2)
1498 :
1499 : ! If not found, skip the loop
1500 45098 : IF (k_start /= 0) THEN
1501 :
1502 207356 : DO k = k_start, k_end
1503 207332 : IF (glob_atm_id(1) < amb_info%raw_torsion_id(1, k)) EXIT
1504 613232 : IF (ANY((glob_atm_id - amb_info%raw_torsion_id(1:4, k)) /= 0)) CYCLE
1505 :
1506 40364 : raw_parm_id = amb_info%raw_torsion_id(5, k)
1507 40364 : imul = torsion_list(j)%torsion_kind%nmul + 1
1508 40364 : CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul)
1509 40364 : CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul)
1510 40364 : CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul)
1511 40364 : torsion_list(j)%torsion_kind%id_type = do_ff_amber
1512 40364 : torsion_list(j)%torsion_kind%k(imul) = amb_info%raw_torsion_k(raw_parm_id)
1513 40364 : torsion_list(j)%torsion_kind%m(imul) = NINT(amb_info%raw_torsion_m(raw_parm_id))
1514 40364 : torsion_list(j)%torsion_kind%phi0(imul) = amb_info%raw_torsion_phi0(raw_parm_id)
1515 40364 : torsion_list(j)%torsion_kind%nmul = imul
1516 207356 : found = .TRUE.
1517 : END DO
1518 :
1519 : END IF
1520 :
1521 : END IF
1522 :
1523 : ! Always have the input param last to overwrite all the other ones
1524 113867 : IF (ASSOCIATED(inp_info%torsion_a)) THEN
1525 192 : DO k = 1, SIZE(inp_info%torsion_a)
1526 : IF ((((inp_info%torsion_a(k)) == (name_atm_a)) .AND. &
1527 : ((inp_info%torsion_b(k)) == (name_atm_b)) .AND. &
1528 : ((inp_info%torsion_c(k)) == (name_atm_c)) .AND. &
1529 166 : ((inp_info%torsion_d(k)) == (name_atm_d))) .OR. &
1530 : (((inp_info%torsion_a(k)) == (name_atm_d)) .AND. &
1531 : ((inp_info%torsion_b(k)) == (name_atm_c)) .AND. &
1532 : ((inp_info%torsion_c(k)) == (name_atm_b)) .AND. &
1533 26 : ((inp_info%torsion_d(k)) == (name_atm_a)))) THEN
1534 38 : imul = torsion_list(j)%torsion_kind%nmul + 1
1535 38 : CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul)
1536 38 : CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul)
1537 38 : CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul)
1538 38 : torsion_list(j)%torsion_kind%id_type = inp_info%torsion_kind(k)
1539 38 : torsion_list(j)%torsion_kind%k(imul) = inp_info%torsion_k(k)
1540 38 : torsion_list(j)%torsion_kind%m(imul) = inp_info%torsion_m(k)
1541 38 : torsion_list(j)%torsion_kind%phi0(imul) = inp_info%torsion_phi0(k)
1542 38 : torsion_list(j)%torsion_kind%nmul = imul
1543 38 : found = .TRUE.
1544 : END IF
1545 : END DO
1546 : END IF
1547 :
1548 113867 : IF (found) THEN
1549 80089 : ldum = cp_to_string(imul)
1550 80089 : IF (iw > 0) THEN
1551 1518 : IF (imul < 1) THEN
1552 : WRITE (UNIT=iw, FMT="(T2,A)") &
1553 0 : "FORCEFIELD| No torsion term found"
1554 1518 : ELSE IF (imul == 1) THEN
1555 : WRITE (UNIT=iw, FMT="(T2,A)") &
1556 : "FORCEFIELD| Found torsion term for the atomic kinds "// &
1557 : TRIM(name_atm_a)//", "//TRIM(name_atm_b)// &
1558 : ", "//TRIM(name_atm_c)// &
1559 1389 : " and "//TRIM(name_atm_d)
1560 : ELSE
1561 : WRITE (UNIT=iw, FMT="(T2,A)") &
1562 : "FORCEFIELD| Found multiple ("//TRIM(ldum)// &
1563 : ") torsion terms for the atomic kinds "// &
1564 : TRIM(name_atm_a)//", "//TRIM(name_atm_b)// &
1565 : ", "//TRIM(name_atm_c)// &
1566 129 : " and "//TRIM(name_atm_d)
1567 : END IF
1568 : END IF
1569 : ELSE
1570 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
1571 : atm2=TRIM(name_atm_b), &
1572 : atm3=TRIM(name_atm_c), &
1573 : atm4=TRIM(name_atm_d), &
1574 : type_name="Torsion", &
1575 33778 : array=Ainfo)
1576 33778 : torsion_list(j)%torsion_kind%id_type = do_ff_undef
1577 33778 : torsion_list(j)%id_type = do_ff_undef
1578 : END IF
1579 :
1580 : ! QM/MM modifications
1581 113867 : IF (only_qm) THEN
1582 1968 : IF (iw > 0) THEN
1583 : WRITE (UNIT=iw, FMT="(T2,A,I0,4(A,I0))") &
1584 0 : "FORCEFIELD| Torsion ", j, " for molecule kind "//TRIM(molecule_name)// &
1585 : TRIM(name_atm_a)// &
1586 : "-"//TRIM(name_atm_b)//"-"//TRIM(name_atm_c)//"-"// &
1587 0 : TRIM(name_atm_d)//" (", torsion_list(j)%a, ", ", &
1588 0 : torsion_list(j)%b, ", ", torsion_list(j)%c, ", ", &
1589 0 : torsion_list(j)%d
1590 : END IF
1591 1968 : torsion_list(j)%torsion_kind%id_type = do_ff_undef
1592 1968 : torsion_list(j)%id_type = do_ff_undef
1593 : END IF
1594 :
1595 : END IF
1596 :
1597 : END DO ! torsion
1598 :
1599 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1600 148287 : torsion_list=torsion_list)
1601 :
1602 : END DO ! molecule kind
1603 :
1604 2643 : CALL timestop(handle2)
1605 :
1606 2643 : END SUBROUTINE force_field_pack_tors
1607 :
1608 : ! **************************************************************************************************
1609 : !> \brief Pack in impropers information needed for the force_field
1610 : !> \param particle_set ...
1611 : !> \param molecule_kind_set ...
1612 : !> \param molecule_set ...
1613 : !> \param Ainfo ...
1614 : !> \param chm_info ...
1615 : !> \param inp_info ...
1616 : !> \param gro_info ...
1617 : !> \param iw ...
1618 : ! **************************************************************************************************
1619 2643 : SUBROUTINE force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, &
1620 : Ainfo, chm_info, inp_info, gro_info, iw)
1621 :
1622 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1623 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1624 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1625 : CHARACTER(LEN=default_string_length), &
1626 : DIMENSION(:), POINTER :: Ainfo
1627 : TYPE(charmm_info_type), POINTER :: chm_info
1628 : TYPE(input_info_type), POINTER :: inp_info
1629 : TYPE(gromos_info_type), POINTER :: gro_info
1630 : INTEGER, INTENT(IN) :: iw
1631 :
1632 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_impr'
1633 :
1634 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c, &
1635 : name_atm_d
1636 : INTEGER :: atm_a, atm_b, atm_c, atm_d, first, &
1637 : handle2, i, itype, j, k, last, natom, &
1638 : nimpr
1639 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1640 : LOGICAL :: found, only_qm
1641 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1642 2643 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list
1643 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1644 : TYPE(molecule_type), POINTER :: molecule
1645 :
1646 2643 : CALL timeset(routineN, handle2)
1647 :
1648 2643 : IF (iw > 0) THEN
1649 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1650 239 : "FORCEFIELD| Checking for improper terms"
1651 : END IF
1652 :
1653 75465 : DO i = 1, SIZE(molecule_kind_set)
1654 :
1655 72822 : molecule_kind => molecule_kind_set(i)
1656 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1657 : molecule_list=molecule_list, &
1658 : natom=natom, &
1659 : nimpr=nimpr, &
1660 72822 : impr_list=impr_list)
1661 :
1662 72822 : molecule => molecule_set(molecule_list(1))
1663 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1664 :
1665 78134 : DO j = 1, nimpr
1666 5312 : atm_a = impr_list(j)%a
1667 5312 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1668 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1669 5312 : name=name_atm_a)
1670 5312 : atm_b = impr_list(j)%b
1671 5312 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1672 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1673 5312 : name=name_atm_b)
1674 5312 : atm_c = impr_list(j)%c
1675 5312 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1676 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1677 5312 : name=name_atm_c)
1678 5312 : atm_d = impr_list(j)%d
1679 5312 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
1680 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1681 5312 : name=name_atm_d)
1682 5312 : found = .FALSE.
1683 5312 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d)
1684 5312 : CALL uppercase(name_atm_a)
1685 5312 : CALL uppercase(name_atm_b)
1686 5312 : CALL uppercase(name_atm_c)
1687 5312 : CALL uppercase(name_atm_d)
1688 :
1689 : ! Loop over params from GROMOS
1690 5312 : IF (ASSOCIATED(gro_info%impr_k)) THEN
1691 0 : k = SIZE(gro_info%impr_k)
1692 0 : itype = impr_list(j)%itype
1693 0 : IF (itype > 0) THEN
1694 0 : impr_list(j)%impr_kind%k = gro_info%impr_k(itype)
1695 0 : impr_list(j)%impr_kind%phi0 = gro_info%impr_phi0(itype)
1696 : ELSE
1697 0 : impr_list(j)%impr_kind%k = gro_info%impr_k(itype)
1698 0 : impr_list(j)%impr_kind%phi0 = gro_info%impr_phi0(itype)
1699 : END IF
1700 0 : found = .TRUE.
1701 0 : impr_list(j)%impr_kind%id_type = gro_info%ff_gromos_type
1702 0 : impr_list(j)%id_type = gro_info%ff_gromos_type
1703 : END IF
1704 :
1705 : ! Loop over params from CHARMM
1706 5312 : IF (ASSOCIATED(chm_info%impr_a)) THEN
1707 171282 : DO k = 1, SIZE(chm_info%impr_a)
1708 : IF ((((chm_info%impr_a(k)) == (name_atm_a)) .AND. &
1709 : ((chm_info%impr_b(k)) == (name_atm_b)) .AND. &
1710 : ((chm_info%impr_c(k)) == (name_atm_c)) .AND. &
1711 168054 : ((chm_info%impr_d(k)) == (name_atm_d))) .OR. &
1712 : (((chm_info%impr_a(k)) == (name_atm_d)) .AND. &
1713 : ((chm_info%impr_b(k)) == (name_atm_c)) .AND. &
1714 : ((chm_info%impr_c(k)) == (name_atm_b)) .AND. &
1715 3228 : ((chm_info%impr_d(k)) == (name_atm_a)))) THEN
1716 1130 : impr_list(j)%impr_kind%id_type = do_ff_charmm
1717 1130 : impr_list(j)%impr_kind%k = chm_info%impr_k(k)
1718 1130 : impr_list(j)%impr_kind%phi0 = chm_info%impr_phi0(k)
1719 : CALL issue_duplications(found, "Impropers", name_atm_a, name_atm_b, &
1720 1130 : name_atm_c, name_atm_d)
1721 1130 : found = .TRUE.
1722 1130 : EXIT
1723 : END IF
1724 : END DO
1725 4358 : IF (.NOT. found) THEN
1726 116678 : DO k = 1, SIZE(chm_info%impr_a)
1727 : IF ((((chm_info%impr_a(k)) == (name_atm_a)) .AND. &
1728 : ((chm_info%impr_b(k)) == ("X")) .AND. &
1729 : ((chm_info%impr_c(k)) == ("X")) .AND. &
1730 115728 : ((chm_info%impr_d(k)) == (name_atm_d))) .OR. &
1731 : (((chm_info%impr_a(k)) == (name_atm_d)) .AND. &
1732 : ((chm_info%impr_b(k)) == ("X")) .AND. &
1733 : ((chm_info%impr_c(k)) == ("X")) .AND. &
1734 950 : ((chm_info%impr_d(k)) == (name_atm_a)))) THEN
1735 2278 : impr_list(j)%impr_kind%id_type = do_ff_charmm
1736 2278 : impr_list(j)%impr_kind%k = chm_info%impr_k(k)
1737 2278 : impr_list(j)%impr_kind%phi0 = chm_info%impr_phi0(k)
1738 : CALL issue_duplications(found, "Impropers", name_atm_a, name_atm_b, &
1739 2278 : name_atm_c, name_atm_d)
1740 2278 : found = .TRUE.
1741 2278 : EXIT
1742 : END IF
1743 : END DO
1744 : END IF
1745 : END IF
1746 :
1747 : ! Loop over params from AMBER not needed since impropers in AMBER
1748 : ! are treated like standard torsions
1749 :
1750 : ! always have the input param last to overwrite all the other ones
1751 5312 : IF (ASSOCIATED(inp_info%impr_a)) THEN
1752 20 : DO k = 1, SIZE(inp_info%impr_a)
1753 : IF (((inp_info%impr_a(k)) == (name_atm_a)) .AND. &
1754 14 : ((inp_info%impr_b(k)) == (name_atm_b)) .AND. &
1755 : ((((inp_info%impr_c(k)) == (name_atm_c)) .AND. &
1756 : ((inp_info%impr_d(k)) == (name_atm_d))) .OR. &
1757 : (((inp_info%impr_c(k)) == (name_atm_d)) .AND. &
1758 6 : ((inp_info%impr_d(k)) == (name_atm_c))))) THEN
1759 8 : impr_list(j)%impr_kind%id_type = inp_info%impr_kind(k)
1760 8 : impr_list(j)%impr_kind%k = inp_info%impr_k(k)
1761 8 : IF (((inp_info%impr_c(k)) == (name_atm_c)) .AND. &
1762 : ((inp_info%impr_d(k)) == (name_atm_d))) THEN
1763 8 : impr_list(j)%impr_kind%phi0 = inp_info%impr_phi0(k)
1764 : ELSE
1765 0 : impr_list(j)%impr_kind%phi0 = -inp_info%impr_phi0(k)
1766 : ! alternative solutions:
1767 : ! - swap impr_list(j)%c with impr_list(j)%d and
1768 : ! name_atom_c with name_atom_d and
1769 : ! atm_c with atm_d
1770 : ! - introduce impr_list(j)%impr_kind%sign. if one, the
1771 : ! sign of phi is not changed in mol_force.f90. if minus
1772 : ! one, the sign of phi is changed in mol_force.f90
1773 : ! similar problems with parameters from charmm pot file
1774 : ! above
1775 : END IF
1776 : CALL issue_duplications(found, "Impropers", name_atm_a, name_atm_b, &
1777 8 : name_atm_c, name_atm_d)
1778 8 : found = .TRUE.
1779 8 : EXIT
1780 : END IF
1781 : END DO
1782 : END IF
1783 :
1784 5312 : IF (.NOT. found) THEN
1785 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
1786 : atm2=TRIM(name_atm_b), &
1787 : atm3=TRIM(name_atm_c), &
1788 : atm4=TRIM(name_atm_d), &
1789 : type_name="Improper", &
1790 1896 : array=Ainfo)
1791 1896 : impr_list(j)%impr_kind%k = 0.0_dp
1792 1896 : impr_list(j)%impr_kind%phi0 = 0.0_dp
1793 1896 : impr_list(j)%impr_kind%id_type = do_ff_undef
1794 1896 : impr_list(j)%id_type = do_ff_undef
1795 : END IF
1796 :
1797 : ! QM/MM modifications
1798 78134 : IF (only_qm) THEN
1799 58 : IF (found) THEN
1800 10 : IF (iw > 0) THEN
1801 : WRITE (UNIT=iw, FMT="(T2,A)") &
1802 : "FORCEFIELD| Found improper term for "//TRIM(name_atm_a)// &
1803 : "-"//TRIM(name_atm_b)//"-"//TRIM(name_atm_c)//"-"// &
1804 0 : TRIM(name_atm_d)
1805 : END IF
1806 : END IF
1807 58 : impr_list(j)%impr_kind%id_type = do_ff_undef
1808 58 : impr_list(j)%id_type = do_ff_undef
1809 : END IF
1810 :
1811 : END DO
1812 :
1813 148287 : CALL set_molecule_kind(molecule_kind=molecule_kind, impr_list=impr_list)
1814 :
1815 : END DO
1816 :
1817 2643 : CALL timestop(handle2)
1818 :
1819 2643 : END SUBROUTINE force_field_pack_impr
1820 :
1821 : ! **************************************************************************************************
1822 : !> \brief Pack in opbend information needed for the force_field.
1823 : !> No loop over params for charmm, amber and gromos since these force
1824 : !> fields have no opbends
1825 : !> \param particle_set ...
1826 : !> \param molecule_kind_set ...
1827 : !> \param molecule_set ...
1828 : !> \param Ainfo ...
1829 : !> \param inp_info ...
1830 : !> \param iw ...
1831 : !> \author Louis Vanduyfhuys
1832 : ! **************************************************************************************************
1833 2643 : SUBROUTINE force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, Ainfo, &
1834 : inp_info, iw)
1835 :
1836 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1837 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1838 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1839 : CHARACTER(LEN=default_string_length), &
1840 : DIMENSION(:), POINTER :: Ainfo
1841 : TYPE(input_info_type), POINTER :: inp_info
1842 : INTEGER, INTENT(IN) :: iw
1843 :
1844 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_opbend'
1845 :
1846 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c, &
1847 : name_atm_d
1848 : INTEGER :: atm_a, atm_b, atm_c, atm_d, first, &
1849 : handle2, i, j, k, last, natom, nopbend
1850 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1851 : LOGICAL :: found, only_qm
1852 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1853 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1854 : TYPE(molecule_type), POINTER :: molecule
1855 2643 : TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list
1856 :
1857 2643 : CALL timeset(routineN, handle2)
1858 :
1859 2643 : IF (iw > 0) THEN
1860 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1861 239 : "FORCEFIELD| Checking for out-of-plane bend terms"
1862 : END IF
1863 :
1864 75465 : DO i = 1, SIZE(molecule_kind_set)
1865 72822 : molecule_kind => molecule_kind_set(i)
1866 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1867 : molecule_list=molecule_list, &
1868 : natom=natom, &
1869 72822 : nopbend=nopbend, opbend_list=opbend_list)
1870 72822 : molecule => molecule_set(molecule_list(1))
1871 :
1872 72822 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1873 78134 : DO j = 1, nopbend
1874 5312 : atm_a = opbend_list(j)%a
1875 5312 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1876 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1877 5312 : name=name_atm_a)
1878 5312 : atm_b = opbend_list(j)%b
1879 5312 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1880 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1881 5312 : name=name_atm_b)
1882 5312 : atm_c = opbend_list(j)%c
1883 5312 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1884 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1885 5312 : name=name_atm_c)
1886 5312 : atm_d = opbend_list(j)%d
1887 5312 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
1888 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1889 5312 : name=name_atm_d)
1890 5312 : found = .FALSE.
1891 5312 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d)
1892 5312 : CALL uppercase(name_atm_a)
1893 5312 : CALL uppercase(name_atm_b)
1894 5312 : CALL uppercase(name_atm_c)
1895 5312 : CALL uppercase(name_atm_d)
1896 :
1897 : ! always have the input param last to overwrite all the other ones
1898 5312 : IF (ASSOCIATED(inp_info%opbend_a)) THEN
1899 2 : DO k = 1, SIZE(inp_info%opbend_a)
1900 : IF (((inp_info%opbend_a(k)) == (name_atm_a)) .AND. &
1901 2 : ((inp_info%opbend_d(k)) == (name_atm_d)) .AND. &
1902 : ((((inp_info%opbend_c(k)) == (name_atm_c)) .AND. &
1903 : ((inp_info%opbend_b(k)) == (name_atm_b))) .OR. &
1904 : (((inp_info%opbend_c(k)) == (name_atm_b)) .AND. &
1905 0 : ((inp_info%opbend_b(k)) == (name_atm_c))))) THEN
1906 2 : opbend_list(j)%opbend_kind%id_type = inp_info%opbend_kind(k)
1907 2 : opbend_list(j)%opbend_kind%k = inp_info%opbend_k(k)
1908 2 : IF (((inp_info%opbend_c(k)) == (name_atm_c)) .AND. &
1909 : ((inp_info%opbend_b(k)) == (name_atm_b))) THEN
1910 2 : opbend_list(j)%opbend_kind%phi0 = inp_info%opbend_phi0(k)
1911 : ELSE
1912 0 : opbend_list(j)%opbend_kind%phi0 = -inp_info%opbend_phi0(k)
1913 : ! alternative solutions:
1914 : ! - swap opbend_list(j)%c with opbend_list(j)%b and
1915 : ! name_atom_c with name_atom_b and
1916 : ! atm_c with atm_b
1917 : ! - introduce opbend_list(j)%opbend_kind%sign. if one, the
1918 : ! sign of phi is not changed in mol_force.f90. if minus
1919 : ! one, the sign of phi is changed in mol_force.f90
1920 :
1921 : END IF
1922 : CALL issue_duplications(found, "Out of plane bend", name_atm_a, name_atm_b, &
1923 2 : name_atm_c, name_atm_d)
1924 2 : found = .TRUE.
1925 2 : EXIT
1926 : END IF
1927 : END DO
1928 : END IF
1929 :
1930 5312 : IF (.NOT. found) THEN
1931 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
1932 : atm2=TRIM(name_atm_b), &
1933 : atm3=TRIM(name_atm_c), &
1934 : atm4=TRIM(name_atm_d), &
1935 : type_name="Out of plane bend", &
1936 5310 : array=Ainfo)
1937 5310 : opbend_list(j)%opbend_kind%k = 0.0_dp
1938 5310 : opbend_list(j)%opbend_kind%phi0 = 0.0_dp
1939 5310 : opbend_list(j)%opbend_kind%id_type = do_ff_undef
1940 5310 : opbend_list(j)%id_type = do_ff_undef
1941 : END IF
1942 : !
1943 : ! QM/MM modifications
1944 : !
1945 78134 : IF (only_qm) THEN
1946 58 : opbend_list(j)%opbend_kind%id_type = do_ff_undef
1947 58 : opbend_list(j)%id_type = do_ff_undef
1948 : END IF
1949 :
1950 : END DO
1951 :
1952 148287 : CALL set_molecule_kind(molecule_kind=molecule_kind, opbend_list=opbend_list)
1953 :
1954 : END DO
1955 :
1956 2643 : CALL timestop(handle2)
1957 :
1958 2643 : END SUBROUTINE force_field_pack_opbend
1959 :
1960 : ! **************************************************************************************************
1961 : !> \brief Set up array of full charges
1962 : !> \param charges ...
1963 : !> \param charges_section ...
1964 : !> \param particle_set ...
1965 : !> \param my_qmmm ...
1966 : !> \param qmmm_env ...
1967 : !> \param inp_info ...
1968 : !> \param iw4 ...
1969 : !> \date 12.2010
1970 : !> \author Teodoro Laino (teodoro.laino@gmail.com)
1971 : ! **************************************************************************************************
1972 8 : SUBROUTINE force_field_pack_charges(charges, charges_section, particle_set, &
1973 : my_qmmm, qmmm_env, inp_info, iw4)
1974 :
1975 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
1976 : TYPE(section_vals_type), POINTER :: charges_section
1977 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1978 : LOGICAL :: my_qmmm
1979 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
1980 : TYPE(input_info_type), POINTER :: inp_info
1981 : INTEGER, INTENT(IN) :: iw4
1982 :
1983 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charges'
1984 :
1985 : CHARACTER(LEN=default_string_length) :: atmname
1986 : INTEGER :: handle, iatom, ilink, j, nval
1987 : LOGICAL :: found_p, is_link_atom, is_ok, &
1988 : only_manybody, only_qm
1989 : REAL(KIND=dp) :: charge, charge_tot, rval, scale_factor
1990 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1991 : TYPE(cp_sll_val_type), POINTER :: list
1992 : TYPE(fist_potential_type), POINTER :: fist_potential
1993 : TYPE(val_type), POINTER :: val
1994 :
1995 8 : CALL timeset(routineN, handle)
1996 :
1997 8 : charge_tot = 0.0_dp
1998 8 : NULLIFY (list)
1999 :
2000 : ! Not implemented for core-shell
2001 8 : IF (ASSOCIATED(inp_info%shell_list)) THEN
2002 0 : CPABORT("Array of charges is not implemented for the core-shell model")
2003 : END IF
2004 :
2005 : ! Allocate array to particle_set size
2006 8 : CPASSERT(.NOT. (ASSOCIATED(charges)))
2007 24 : ALLOCATE (charges(SIZE(particle_set)))
2008 :
2009 : ! Fill with input values
2010 8 : CALL section_vals_val_get(charges_section, "_DEFAULT_KEYWORD_", n_rep_val=nval)
2011 8 : CPASSERT(nval == SIZE(charges))
2012 8 : CALL section_vals_list_get(charges_section, "_DEFAULT_KEYWORD_", list=list)
2013 44 : DO iatom = 1, nval
2014 : ! we use only the first default_string_length characters of each line
2015 36 : is_ok = cp_sll_val_next(list, val)
2016 36 : CALL val_get(val, r_val=rval)
2017 : ! assign values
2018 36 : charges(iatom) = rval
2019 :
2020 : ! Perform a post-processing
2021 36 : atomic_kind => particle_set(iatom)%atomic_kind
2022 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2023 : fist_potential=fist_potential, &
2024 36 : name=atmname)
2025 36 : CALL get_potential(potential=fist_potential, qeff=charge)
2026 :
2027 36 : only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom)
2028 36 : CALL uppercase(atmname)
2029 36 : IF (charge /= -HUGE(0.0_dp)) &
2030 : CALL cp_warn(__LOCATION__, &
2031 : "The charge for atom index ("//cp_to_string(iatom)//") and atom name ("// &
2032 : TRIM(atmname)//") was already defined. The charge associated to this kind"// &
2033 0 : " will be set to an uninitialized value and only the atom specific charge will be used! ")
2034 36 : charge = -HUGE(0.0_dp)
2035 :
2036 : ! Check if the potential really requires the charge definition..
2037 36 : IF (ASSOCIATED(inp_info%nonbonded)) THEN
2038 18 : IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN
2039 : ! Let's find the nonbonded potential where this atom is involved
2040 18 : only_manybody = .TRUE.
2041 18 : found_p = .FALSE.
2042 30 : DO j = 1, SIZE(inp_info%nonbonded%pot)
2043 30 : IF (atmname == inp_info%nonbonded%pot(j)%pot%at1 .OR. &
2044 0 : atmname == inp_info%nonbonded%pot(j)%pot%at2) THEN
2045 18 : SELECT CASE (inp_info%nonbonded%pot(j)%pot%type(1))
2046 : CASE (ea_type, tersoff_type, siepmann_type)
2047 : ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential
2048 : ! Do nothing..
2049 : CASE DEFAULT
2050 : only_manybody = .FALSE.
2051 18 : EXIT
2052 : END SELECT
2053 : found_p = .TRUE.
2054 : END IF
2055 : END DO
2056 18 : IF (only_manybody .AND. found_p) THEN
2057 0 : charges(iatom) = 0.0_dp
2058 : END IF
2059 : END IF
2060 : END IF
2061 :
2062 : ! QM/MM modifications
2063 80 : IF (only_qm .AND. my_qmmm) THEN
2064 6 : IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN
2065 6 : scale_factor = 0.0_dp
2066 6 : IF (is_link_atom) THEN
2067 : ! Find the scaling factor...
2068 0 : DO ilink = 1, SIZE(qmmm_env%mm_link_atoms)
2069 0 : IF (iatom == qmmm_env%mm_link_atoms(ilink)) EXIT
2070 : END DO
2071 0 : CPASSERT(ilink <= SIZE(qmmm_env%mm_link_atoms))
2072 0 : scale_factor = qmmm_env%fist_scale_charge_link(ilink)
2073 : END IF
2074 6 : charges(iatom) = charges(iatom)*scale_factor
2075 : END IF
2076 : END IF
2077 : END DO
2078 :
2079 : ! Sum up total charges for IO
2080 44 : charge_tot = SUM(charges)
2081 :
2082 : ! Print total charge of the system
2083 8 : IF (iw4 > 0) THEN
2084 : WRITE (UNIT=iw4, FMT="(/,T2,A,T61,F20.10)") &
2085 4 : "FORCEFIELD| Total charge of the classical system: ", charge_tot
2086 : END IF
2087 :
2088 8 : CALL timestop(handle)
2089 :
2090 16 : END SUBROUTINE force_field_pack_charges
2091 :
2092 : ! **************************************************************************************************
2093 : !> \brief Set up atomic_kind_set()%fist_potential%[qeff]
2094 : !> and shell potential parameters
2095 : !> \param atomic_kind_set ...
2096 : !> \param qmmm_env ...
2097 : !> \param fatal ...
2098 : !> \param iw ...
2099 : !> \param iw4 ...
2100 : !> \param Ainfo ...
2101 : !> \param my_qmmm ...
2102 : !> \param inp_info ...
2103 : ! **************************************************************************************************
2104 2635 : SUBROUTINE force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4, &
2105 : Ainfo, my_qmmm, inp_info)
2106 :
2107 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2108 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
2109 : LOGICAL, INTENT(INOUT) :: fatal
2110 : INTEGER, INTENT(IN) :: iw, iw4
2111 : CHARACTER(LEN=default_string_length), &
2112 : DIMENSION(:), POINTER :: Ainfo
2113 : LOGICAL, INTENT(IN) :: my_qmmm
2114 : TYPE(input_info_type), POINTER :: inp_info
2115 :
2116 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charge'
2117 :
2118 : CHARACTER(LEN=default_string_length) :: atmname
2119 : INTEGER :: handle, i, ilink, j
2120 2635 : INTEGER, DIMENSION(:), POINTER :: my_atom_list
2121 : LOGICAL :: found, found_p, is_link_atom, is_shell, &
2122 : only_manybody, only_qm
2123 : REAL(KIND=dp) :: charge, charge_tot, cs_charge, &
2124 : scale_factor
2125 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2126 : TYPE(fist_potential_type), POINTER :: fist_potential
2127 :
2128 2635 : CALL timeset(routineN, handle)
2129 :
2130 2635 : charge_tot = 0.0_dp
2131 :
2132 13894 : DO i = 1, SIZE(atomic_kind_set)
2133 11259 : atomic_kind => atomic_kind_set(i)
2134 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2135 : fist_potential=fist_potential, &
2136 : atom_list=my_atom_list, &
2137 11259 : name=atmname)
2138 11259 : CALL get_potential(potential=fist_potential, qeff=charge)
2139 :
2140 11259 : is_shell = .FALSE.
2141 11259 : found = .FALSE.
2142 11259 : only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom)
2143 11259 : CALL uppercase(atmname)
2144 11259 : IF (charge /= -HUGE(0.0_dp)) found = .TRUE.
2145 :
2146 : ! Always have the input param last to overwrite all the other ones
2147 11259 : IF (ASSOCIATED(inp_info%charge_atm)) THEN
2148 5581 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") ""
2149 27332 : DO j = 1, SIZE(inp_info%charge_atm)
2150 : IF (debug_this_module) THEN
2151 : IF (iw > 0) THEN
2152 : WRITE (UNIT=iw, FMT="(T2,A)") &
2153 : "Checking charges for the atomic kinds "// &
2154 : TRIM(inp_info%charge_atm(j))//" and "//TRIM(atmname)
2155 : END IF
2156 : END IF
2157 27332 : IF ((inp_info%charge_atm(j)) == atmname) THEN
2158 5483 : charge = inp_info%charge(j)
2159 5483 : CALL issue_duplications(found, "Charge", atmname)
2160 5483 : found = .TRUE.
2161 : END IF
2162 : END DO
2163 : END IF
2164 : ! Check if the ATOM type has a core-shell associated.. In this case
2165 : ! print a warning: the CHARGE will not be used if defined.. or we can
2166 : ! even skip its definition..
2167 11259 : IF (ASSOCIATED(inp_info%shell_list)) THEN
2168 1414 : DO j = 1, SIZE(inp_info%shell_list)
2169 1414 : IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN
2170 448 : is_shell = .TRUE.
2171 : cs_charge = inp_info%shell_list(j)%shell%charge_core + &
2172 448 : inp_info%shell_list(j)%shell%charge_shell
2173 448 : charge = 0.0_dp
2174 448 : IF (found) THEN
2175 : IF (found) &
2176 : CALL cp_warn(__LOCATION__, &
2177 : "CORE-SHELL model defined for KIND ("//TRIM(atmname)//")"// &
2178 204 : " ignoring charge definition! ")
2179 : ELSE
2180 244 : found = .TRUE.
2181 : END IF
2182 : END IF
2183 : END DO
2184 : END IF
2185 : ! Check if the potential really requires the charge definition..
2186 11259 : IF (ASSOCIATED(inp_info%nonbonded)) THEN
2187 4347 : IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN
2188 : ! Let's find the nonbonded potential where this atom is involved
2189 4347 : only_manybody = .TRUE.
2190 4347 : found_p = .FALSE.
2191 7881 : DO j = 1, SIZE(inp_info%nonbonded%pot)
2192 7678 : IF (atmname == inp_info%nonbonded%pot(j)%pot%at1 .OR. &
2193 203 : atmname == inp_info%nonbonded%pot(j)%pot%at2) THEN
2194 4340 : SELECT CASE (inp_info%nonbonded%pot(j)%pot%type(1))
2195 : CASE (ea_type, tersoff_type, siepmann_type, quip_type, nequip_type, &
2196 : allegro_type, deepmd_type, ace_type)
2197 : ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential
2198 : ! Do nothing..
2199 : CASE DEFAULT
2200 : only_manybody = .FALSE.
2201 4340 : EXIT
2202 : END SELECT
2203 : found_p = .TRUE.
2204 : END IF
2205 : END DO
2206 4347 : IF (only_manybody .AND. found_p) THEN
2207 156 : charge = 0.0_dp
2208 156 : found = .TRUE.
2209 : END IF
2210 : END IF
2211 : END IF
2212 11259 : IF (.NOT. found) THEN
2213 : ! Set the charge to zero anyway in case the user decides to ignore
2214 : ! missing critical parameters.
2215 12 : charge = 0.0_dp
2216 : CALL store_FF_missing_par(atm1=TRIM(atmname), &
2217 : fatal=fatal, &
2218 : type_name="Charge", &
2219 12 : array=Ainfo)
2220 : END IF
2221 : !
2222 : ! QM/MM modifications
2223 : !
2224 11259 : IF (only_qm .AND. my_qmmm) THEN
2225 1286 : IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN
2226 1076 : scale_factor = 0.0_dp
2227 1076 : IF (is_link_atom) THEN
2228 : !
2229 : ! Find the scaling factor...
2230 : !
2231 386 : DO ilink = 1, SIZE(qmmm_env%mm_link_atoms)
2232 658 : IF (ANY(my_atom_list == qmmm_env%mm_link_atoms(ilink))) EXIT
2233 : END DO
2234 114 : CPASSERT(ilink <= SIZE(qmmm_env%mm_link_atoms))
2235 114 : scale_factor = qmmm_env%fist_scale_charge_link(ilink)
2236 : END IF
2237 1076 : charge = charge*scale_factor
2238 : END IF
2239 : END IF
2240 :
2241 11259 : CALL set_potential(potential=fist_potential, qeff=charge)
2242 : ! Sum up total charges for IO
2243 13894 : IF (found) THEN
2244 11247 : IF (is_shell) THEN
2245 448 : charge_tot = charge_tot + atomic_kind%natom*cs_charge
2246 : ELSE
2247 10799 : charge_tot = charge_tot + atomic_kind%natom*charge
2248 : END IF
2249 : END IF
2250 : END DO
2251 :
2252 : ! Print total charge of the system
2253 2635 : IF (iw4 > 0) THEN
2254 : WRITE (UNIT=iw4, FMT="(/,T2,A,T61,F20.10)") &
2255 1301 : "FORCEFIELD| Total charge of the classical system: ", charge_tot
2256 : END IF
2257 :
2258 2635 : CALL timestop(handle)
2259 :
2260 2635 : END SUBROUTINE force_field_pack_charge
2261 :
2262 : ! **************************************************************************************************
2263 : !> \brief Set up the radius of the electrostatic multipole in Fist
2264 : !> \param atomic_kind_set ...
2265 : !> \param iw ...
2266 : !> \param subsys_section ...
2267 : !> \author Toon.Verstraelen@gmail.com
2268 : ! **************************************************************************************************
2269 5286 : SUBROUTINE force_field_pack_radius(atomic_kind_set, iw, subsys_section)
2270 :
2271 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2272 : INTEGER, INTENT(IN) :: iw
2273 : TYPE(section_vals_type), POINTER :: subsys_section
2274 :
2275 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_radius'
2276 :
2277 : CHARACTER(LEN=default_string_length) :: inp_kind_name, kind_name
2278 : INTEGER :: handle, i, i_rep, n_rep
2279 : LOGICAL :: found
2280 : REAL(KIND=dp) :: mm_radius
2281 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2282 : TYPE(fist_potential_type), POINTER :: fist_potential
2283 : TYPE(section_vals_type), POINTER :: kind_section
2284 :
2285 2643 : CALL timeset(routineN, handle)
2286 :
2287 2643 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
2288 2643 : CALL section_vals_get(kind_section, n_repetition=n_rep)
2289 :
2290 13922 : DO i = 1, SIZE(atomic_kind_set)
2291 11279 : atomic_kind => atomic_kind_set(i)
2292 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2293 11279 : fist_potential=fist_potential, name=kind_name)
2294 11279 : CALL uppercase(kind_name)
2295 11279 : found = .FALSE.
2296 :
2297 : ! Try to find a matching KIND section in the SUBSYS section and read the
2298 : ! MM_RADIUS field if it is present. In case the kind section is never
2299 : ! encountered, the mm_radius remains zero.
2300 11279 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") ""
2301 11279 : mm_radius = 0.0_dp
2302 39568 : DO i_rep = 1, n_rep
2303 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
2304 28289 : c_val=inp_kind_name, i_rep_section=i_rep)
2305 28289 : CALL uppercase(inp_kind_name)
2306 28289 : IF (iw > 0) THEN
2307 : WRITE (UNIT=iw, FMT="(T2,A)") &
2308 : "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// &
2309 905 : " and "//TRIM(inp_kind_name)//" for MM_RADIUS"
2310 : END IF
2311 39568 : IF (TRIM(kind_name) == TRIM(inp_kind_name)) THEN
2312 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
2313 1839 : keyword_name="MM_RADIUS", r_val=mm_radius)
2314 1839 : CALL issue_duplications(found, "MM_RADIUS", kind_name)
2315 1839 : found = .TRUE.
2316 : END IF
2317 : END DO
2318 13922 : CALL set_potential(potential=fist_potential, mm_radius=mm_radius)
2319 : END DO
2320 :
2321 2643 : CALL timestop(handle)
2322 :
2323 2643 : END SUBROUTINE force_field_pack_radius
2324 :
2325 : ! **************************************************************************************************
2326 : !> \brief Set up the polarizable FF parameters
2327 : !> \param atomic_kind_set ...
2328 : !> \param iw ...
2329 : !> \param inp_info ...
2330 : !> \author Toon.Verstraelen@gmail.com
2331 : ! **************************************************************************************************
2332 2643 : SUBROUTINE force_field_pack_pol(atomic_kind_set, iw, inp_info)
2333 :
2334 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2335 : INTEGER, INTENT(IN) :: iw
2336 : TYPE(input_info_type), POINTER :: inp_info
2337 :
2338 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_pol'
2339 :
2340 : CHARACTER(LEN=default_string_length) :: kind_name
2341 : INTEGER :: handle, i, j
2342 : LOGICAL :: found
2343 : REAL(KIND=dp) :: apol, cpol
2344 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2345 : TYPE(fist_potential_type), POINTER :: fist_potential
2346 :
2347 2643 : CALL timeset(routineN, handle)
2348 :
2349 2643 : IF (iw > 0) THEN
2350 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2351 239 : "FORCEFIELD| Checking for polarisable forcefield terms"
2352 : END IF
2353 :
2354 13922 : DO i = 1, SIZE(atomic_kind_set)
2355 11279 : atomic_kind => atomic_kind_set(i)
2356 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2357 : fist_potential=fist_potential, &
2358 11279 : name=kind_name)
2359 11279 : CALL get_potential(potential=fist_potential, apol=apol, cpol=cpol)
2360 11279 : CALL uppercase(kind_name)
2361 11279 : found = .FALSE.
2362 :
2363 11279 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") ""
2364 : ! Always have the input param last to overwrite all the other ones
2365 11279 : IF (ASSOCIATED(inp_info%apol_atm)) THEN
2366 292 : DO j = 1, SIZE(inp_info%apol_atm)
2367 200 : IF (iw > 0) THEN
2368 : WRITE (UNIT=iw, FMT="(T2,A)") &
2369 : "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// &
2370 0 : " and "//TRIM(inp_info%apol_atm(j))//" for APOL"
2371 : END IF
2372 292 : IF ((inp_info%apol_atm(j)) == kind_name) THEN
2373 64 : apol = inp_info%apol(j)
2374 64 : CALL issue_duplications(found, "APOL", kind_name)
2375 64 : found = .TRUE.
2376 : END IF
2377 : END DO
2378 : END IF
2379 :
2380 11279 : IF (ASSOCIATED(inp_info%cpol_atm)) THEN
2381 0 : DO j = 1, SIZE(inp_info%cpol_atm)
2382 0 : IF (iw > 0) THEN
2383 : WRITE (UNIT=iw, FMT="(T2,A)") &
2384 : "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// &
2385 0 : " and "//TRIM(inp_info%cpol_atm(j))//" for CPOL"
2386 : END IF
2387 0 : IF ((inp_info%cpol_atm(j)) == kind_name) THEN
2388 0 : cpol = inp_info%cpol(j)
2389 0 : CALL issue_duplications(found, "CPOL", kind_name)
2390 0 : found = .TRUE.
2391 : END IF
2392 : END DO
2393 : END IF
2394 :
2395 13922 : CALL set_potential(potential=fist_potential, apol=apol, cpol=cpol)
2396 :
2397 : END DO
2398 :
2399 2643 : CALL timestop(handle)
2400 :
2401 2643 : END SUBROUTINE force_field_pack_pol
2402 :
2403 : ! **************************************************************************************************
2404 : !> \brief Set up damping parameters
2405 : !> \param atomic_kind_set ...
2406 : !> \param iw ...
2407 : !> \param inp_info ...
2408 : ! **************************************************************************************************
2409 2643 : SUBROUTINE force_field_pack_damp(atomic_kind_set, iw, inp_info)
2410 :
2411 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2412 : INTEGER :: iw
2413 : TYPE(input_info_type), POINTER :: inp_info
2414 :
2415 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_damp'
2416 :
2417 : CHARACTER(len=default_string_length) :: atm_name1, atm_name2, my_atm_name1, &
2418 : my_atm_name2
2419 : INTEGER :: handle2, i, j, k, nkinds
2420 : LOGICAL :: found
2421 : TYPE(atomic_kind_type), POINTER :: atomic_kind, atomic_kind2
2422 : TYPE(damping_p_type), POINTER :: damping
2423 :
2424 2643 : CALL timeset(routineN, handle2)
2425 :
2426 2643 : IF (iw > 0) THEN
2427 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2428 239 : "FORCEFIELD| Checking for damping terms"
2429 : END IF
2430 :
2431 2643 : NULLIFY (damping)
2432 2643 : nkinds = SIZE(atomic_kind_set)
2433 :
2434 13922 : DO j = 1, SIZE(atomic_kind_set)
2435 :
2436 11279 : atomic_kind => atomic_kind_set(j)
2437 :
2438 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2439 11279 : name=atm_name1)
2440 11279 : CALL uppercase(atm_name1)
2441 :
2442 11279 : IF (ASSOCIATED(inp_info%damping_list)) THEN
2443 50 : DO i = 1, SIZE(inp_info%damping_list)
2444 28 : my_atm_name1 = inp_info%damping_list(i)%atm_name1
2445 28 : my_atm_name2 = inp_info%damping_list(i)%atm_name2
2446 : IF (debug_this_module) THEN
2447 : IF (iw > 0) THEN
2448 : WRITE (UNIT=iw, FMT="(T2,A)") &
2449 : "FORCEFIELD| Check damping for the atomic kinds "// &
2450 : TRIM(my_atm_name1)//" and "//TRIM(atm_name1)
2451 : END IF
2452 : END IF
2453 50 : IF (my_atm_name1 == atm_name1) THEN
2454 12 : IF (.NOT. ASSOCIATED(damping)) THEN
2455 10 : CALL damping_p_create(damping, nkinds)
2456 : END IF
2457 12 : found = .FALSE.
2458 40 : DO k = 1, SIZE(atomic_kind_set)
2459 28 : atomic_kind2 => atomic_kind_set(k)
2460 : CALL get_atomic_kind(atomic_kind=atomic_kind2, &
2461 28 : name=atm_name2)
2462 28 : CALL uppercase(atm_name2)
2463 40 : IF (my_atm_name2 == atm_name2) THEN
2464 12 : IF (damping%damp(k)%bij /= HUGE(0.0_dp)) found = .TRUE.
2465 12 : CALL issue_duplications(found, "Damping", atm_name1)
2466 12 : found = .TRUE.
2467 24 : SELECT CASE (TRIM(inp_info%damping_list(i)%dtype))
2468 : CASE ('TANG-TOENNIES')
2469 12 : damping%damp(k)%itype = tang_toennies
2470 : CASE DEFAULT
2471 24 : CPABORT("Unknown damping type.")
2472 : END SELECT
2473 12 : damping%damp(k)%order = inp_info%damping_list(i)%order
2474 12 : damping%damp(k)%bij = inp_info%damping_list(i)%bij
2475 12 : damping%damp(k)%cij = inp_info%damping_list(i)%cij
2476 : END IF
2477 : END DO
2478 12 : IF (.NOT. found) THEN
2479 : CALL cp_warn(__LOCATION__, &
2480 : "Atom "//TRIM(my_atm_name2)// &
2481 : " in damping parameters for atom "//TRIM(my_atm_name1)// &
2482 0 : " not found.")
2483 : END IF
2484 : END IF
2485 : END DO
2486 : END IF
2487 :
2488 11279 : CALL set_atomic_kind(atomic_kind=atomic_kind, damping=damping)
2489 :
2490 13922 : NULLIFY (damping)
2491 :
2492 : END DO
2493 :
2494 2643 : CALL timestop(handle2)
2495 :
2496 2643 : END SUBROUTINE force_field_pack_damp
2497 :
2498 : ! **************************************************************************************************
2499 : !> \brief Set up shell potential parameters
2500 : !> \param particle_set ...
2501 : !> \param atomic_kind_set ...
2502 : !> \param molecule_kind_set ...
2503 : !> \param molecule_set ...
2504 : !> \param root_section ...
2505 : !> \param subsys_section ...
2506 : !> \param shell_particle_set ...
2507 : !> \param core_particle_set ...
2508 : !> \param cell ...
2509 : !> \param iw ...
2510 : !> \param inp_info ...
2511 : ! **************************************************************************************************
2512 13215 : SUBROUTINE force_field_pack_shell(particle_set, atomic_kind_set, &
2513 : molecule_kind_set, molecule_set, root_section, subsys_section, &
2514 : shell_particle_set, core_particle_set, cell, iw, inp_info)
2515 :
2516 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2517 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2518 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2519 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
2520 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
2521 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_particle_set, core_particle_set
2522 : TYPE(cell_type), POINTER :: cell
2523 : INTEGER :: iw
2524 : TYPE(input_info_type), POINTER :: inp_info
2525 :
2526 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_shell'
2527 :
2528 : CHARACTER(LEN=default_string_length) :: atmname
2529 : INTEGER :: counter, first, first_shell, handle2, i, &
2530 : j, last, last_shell, n, natom, nmol, &
2531 : nshell_tot
2532 2643 : INTEGER, DIMENSION(:), POINTER :: molecule_list, shell_list_tmp
2533 : LOGICAL :: core_coord_read, found_shell, is_a_shell, is_link_atom, null_massfrac, only_qm, &
2534 : save_mem, shell_adiabatic, shell_coord_read
2535 : REAL(KIND=dp) :: atmmass
2536 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2537 : TYPE(molecule_kind_type), POINTER :: molecule_kind
2538 : TYPE(molecule_type), POINTER :: molecule
2539 : TYPE(section_vals_type), POINTER :: global_section
2540 : TYPE(shell_kind_type), POINTER :: shell
2541 2643 : TYPE(shell_type), DIMENSION(:), POINTER :: shell_list
2542 :
2543 2643 : CALL timeset(routineN, handle2)
2544 :
2545 2643 : nshell_tot = 0
2546 2643 : n = 0
2547 2643 : first_shell = 1
2548 2643 : null_massfrac = .FALSE.
2549 2643 : core_coord_read = .FALSE.
2550 2643 : shell_coord_read = .FALSE.
2551 :
2552 2643 : NULLIFY (global_section)
2553 2643 : global_section => section_vals_get_subs_vals(root_section, "GLOBAL")
2554 2643 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
2555 :
2556 2643 : IF (iw > 0) THEN
2557 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2558 239 : "FORCEFIELD| Checking for core-shell terms"
2559 : END IF
2560 :
2561 13922 : DO i = 1, SIZE(atomic_kind_set)
2562 11279 : atomic_kind => atomic_kind_set(i)
2563 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2564 11279 : name=atmname)
2565 :
2566 11279 : found_shell = .FALSE.
2567 11279 : only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom)
2568 11279 : CALL uppercase(atmname)
2569 :
2570 : ! The shell potential can be defined only from input
2571 13922 : IF (ASSOCIATED(inp_info%shell_list)) THEN
2572 1414 : DO j = 1, SIZE(inp_info%shell_list)
2573 : IF (debug_this_module) THEN
2574 : IF (iw > 0) THEN
2575 : WRITE (UNIT=iw, FMT="(T2,A)") &
2576 : "Checking shells for the atomic kinds "// &
2577 : TRIM(inp_info%shell_list(j)%atm_name)//" and "//TRIM(atmname)
2578 : END IF
2579 : END IF
2580 1414 : IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN
2581 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2582 448 : shell=shell, mass=atmmass, natom=natom)
2583 448 : IF (.NOT. ASSOCIATED(shell)) ALLOCATE (shell)
2584 448 : nshell_tot = nshell_tot + natom
2585 448 : shell%charge_core = inp_info%shell_list(j)%shell%charge_core
2586 448 : shell%charge_shell = inp_info%shell_list(j)%shell%charge_shell
2587 448 : shell%massfrac = inp_info%shell_list(j)%shell%massfrac
2588 448 : IF (shell%massfrac < EPSILON(1.0_dp)) null_massfrac = .TRUE.
2589 448 : shell%k2_spring = inp_info%shell_list(j)%shell%k2_spring
2590 448 : shell%k4_spring = inp_info%shell_list(j)%shell%k4_spring
2591 448 : shell%max_dist = inp_info%shell_list(j)%shell%max_dist
2592 448 : shell%shell_cutoff = inp_info%shell_list(j)%shell%shell_cutoff
2593 448 : shell%mass_shell = shell%massfrac*atmmass
2594 448 : shell%mass_core = atmmass - shell%mass_shell
2595 448 : CALL issue_duplications(found_shell, "Shell", atmname)
2596 448 : found_shell = .TRUE.
2597 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
2598 448 : shell=shell, shell_active=.TRUE.)
2599 : END IF
2600 : END DO ! shell kind
2601 : END IF ! associated shell_list
2602 : END DO ! atomic kind
2603 :
2604 2643 : IF (iw > 0) THEN
2605 : WRITE (UNIT=iw, FMT="(/,T2,A,T61,I20)") &
2606 239 : "FORCEFIELD| Total number of particles with a shell:", nshell_tot
2607 : END IF
2608 : ! If shell-model is present: Create particle_set of shells (coord. vel. force)
2609 2643 : NULLIFY (shell_particle_set)
2610 2643 : NULLIFY (core_particle_set)
2611 2643 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, shell_adiabatic=shell_adiabatic)
2612 2643 : IF (nshell_tot > 0) THEN
2613 256 : IF (shell_adiabatic .AND. null_massfrac) THEN
2614 0 : CPABORT("Shell-model adiabatic: at least one shell_kind has mass zero")
2615 : END IF
2616 256 : CALL allocate_particle_set(shell_particle_set, nshell_tot)
2617 256 : CALL allocate_particle_set(core_particle_set, nshell_tot)
2618 256 : counter = 0
2619 : ! Initialise the shell (and core) coordinates with the particle (atomic) coordinates,
2620 : ! count the shell and set pointers
2621 28934 : DO i = 1, SIZE(particle_set)
2622 28678 : NULLIFY (atomic_kind)
2623 28678 : NULLIFY (shell)
2624 28678 : atomic_kind => particle_set(i)%atomic_kind
2625 28678 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_a_shell)
2626 28934 : IF (is_a_shell) THEN
2627 28070 : counter = counter + 1
2628 28070 : particle_set(i)%shell_index = counter
2629 28070 : shell_particle_set(counter)%shell_index = counter
2630 28070 : shell_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind
2631 196490 : shell_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3)
2632 28070 : shell_particle_set(counter)%atom_index = i
2633 28070 : core_particle_set(counter)%shell_index = counter
2634 28070 : core_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind
2635 196490 : core_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3)
2636 28070 : core_particle_set(counter)%atom_index = i
2637 : ELSE
2638 608 : particle_set(i)%shell_index = 0
2639 : END IF
2640 : END DO
2641 256 : CPASSERT(counter == nshell_tot)
2642 : END IF
2643 :
2644 : ! Read the shell (and core) coordinates from the restart file, if available
2645 : CALL read_binary_cs_coordinates("SHELL", shell_particle_set, root_section, &
2646 2643 : subsys_section, shell_coord_read)
2647 : CALL read_binary_cs_coordinates("CORE", core_particle_set, root_section, &
2648 2643 : subsys_section, core_coord_read)
2649 :
2650 2643 : IF (nshell_tot > 0) THEN
2651 : ! Read the shell (and core) coordinates from the input, if no coordinates were found
2652 : ! in the restart file
2653 256 : IF (shell_adiabatic) THEN
2654 256 : IF (.NOT. (core_coord_read .AND. shell_coord_read)) THEN
2655 : CALL read_shell_coord_input(particle_set, shell_particle_set, cell, &
2656 : subsys_section, core_particle_set, &
2657 240 : save_mem=save_mem)
2658 : END IF
2659 : ELSE
2660 0 : IF (.NOT. shell_coord_read) THEN
2661 : CALL read_shell_coord_input(particle_set, shell_particle_set, cell, &
2662 0 : subsys_section, save_mem=save_mem)
2663 : END IF
2664 : END IF
2665 : ! Determine the number of shells per molecule kind
2666 256 : n = 0
2667 11548 : DO i = 1, SIZE(molecule_kind_set)
2668 11292 : molecule_kind => molecule_kind_set(i)
2669 : CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, &
2670 11292 : natom=natom, nmolecule=nmol)
2671 11292 : molecule => molecule_set(molecule_list(1))
2672 11292 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
2673 33876 : ALLOCATE (shell_list_tmp(natom))
2674 11292 : counter = 0
2675 23562 : DO j = first, last
2676 12270 : atomic_kind => particle_set(j)%atomic_kind
2677 12270 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_a_shell)
2678 23562 : IF (is_a_shell) THEN
2679 11804 : counter = counter + 1
2680 11804 : shell_list_tmp(counter) = j - first + 1
2681 11804 : first_shell = MIN(first_shell, MAX(1, particle_set(j)%shell_index))
2682 : END IF
2683 : END DO ! j atom in molecule_kind i, molecule 1 of the molecule_list
2684 11292 : IF (counter /= 0) THEN
2685 : ! Setup of fist_shell and last_shell for all molecules..
2686 29288 : DO j = 1, SIZE(molecule_list)
2687 18414 : last_shell = first_shell + counter - 1
2688 18414 : molecule => molecule_set(molecule_list(j))
2689 18414 : molecule%first_shell = first_shell
2690 18414 : molecule%last_shell = last_shell
2691 29288 : first_shell = last_shell + 1
2692 : END DO
2693 : ! Setup of shell_list
2694 10874 : CALL get_molecule_kind(molecule_kind=molecule_kind, shell_list=shell_list)
2695 10874 : IF (ASSOCIATED(shell_list)) THEN
2696 0 : DEALLOCATE (shell_list)
2697 : END IF
2698 44426 : ALLOCATE (shell_list(counter))
2699 22678 : DO j = 1, counter
2700 11804 : shell_list(j)%a = shell_list_tmp(j)
2701 11804 : atomic_kind => particle_set(shell_list_tmp(j) + first - 1)%atomic_kind
2702 11804 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname, shell=shell)
2703 11804 : CALL uppercase(atmname)
2704 11804 : shell_list(j)%name = atmname
2705 22678 : shell_list(j)%shell_kind => shell
2706 : END DO
2707 10874 : CALL set_molecule_kind(molecule_kind=molecule_kind, nshell=counter, shell_list=shell_list)
2708 : END IF
2709 11292 : DEALLOCATE (shell_list_tmp)
2710 22840 : n = n + nmol*counter
2711 : END DO ! i molecule kind
2712 : END IF
2713 :
2714 2643 : CPASSERT(first_shell - 1 == nshell_tot)
2715 2643 : CPASSERT(n == nshell_tot)
2716 :
2717 2643 : CALL timestop(handle2)
2718 :
2719 2643 : END SUBROUTINE force_field_pack_shell
2720 :
2721 : ! **************************************************************************************************
2722 : !> \brief Assign input and potential info to potparm_nonbond14
2723 : !> \param atomic_kind_set ...
2724 : !> \param ff_type ...
2725 : !> \param qmmm_env ...
2726 : !> \param iw ...
2727 : !> \param Ainfo ...
2728 : !> \param chm_info ...
2729 : !> \param inp_info ...
2730 : !> \param gro_info ...
2731 : !> \param amb_info ...
2732 : !> \param potparm_nonbond14 ...
2733 : !> \param ewald_env ...
2734 : ! **************************************************************************************************
2735 2627 : SUBROUTINE force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, &
2736 : Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env)
2737 :
2738 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2739 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
2740 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
2741 : INTEGER :: iw
2742 : CHARACTER(LEN=default_string_length), &
2743 : DIMENSION(:), POINTER :: Ainfo
2744 : TYPE(charmm_info_type), POINTER :: chm_info
2745 : TYPE(input_info_type), POINTER :: inp_info
2746 : TYPE(gromos_info_type), POINTER :: gro_info
2747 : TYPE(amber_info_type), POINTER :: amb_info
2748 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond14
2749 : TYPE(ewald_environment_type), POINTER :: ewald_env
2750 :
2751 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond14'
2752 :
2753 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
2754 : name_atm_b, name_atm_b_local
2755 : INTEGER :: handle2, i, ii, j, jj, k, match_names
2756 : LOGICAL :: found, found_a, found_b, only_qm, &
2757 : use_qmmm_ff
2758 : REAL(KIND=dp) :: epsilon0, epsilon_a, epsilon_b, &
2759 : ewald_rcut, rmin, rmin2_a, rmin2_b
2760 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2761 : TYPE(pair_potential_single_type), POINTER :: pot
2762 :
2763 2627 : CALL timeset(routineN, handle2)
2764 :
2765 2627 : use_qmmm_ff = qmmm_env%use_qmmm_ff
2766 2627 : NULLIFY (pot)
2767 :
2768 2627 : IF (iw > 0) THEN
2769 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2770 238 : "FORCEFIELD| Checking for nonbonded14 terms"
2771 : END IF
2772 :
2773 2627 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
2774 2627 : CALL pair_potential_pp_create(potparm_nonbond14, SIZE(atomic_kind_set))
2775 :
2776 13832 : DO i = 1, SIZE(atomic_kind_set)
2777 11205 : atomic_kind => atomic_kind_set(i)
2778 11205 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
2779 271546 : DO j = i, SIZE(atomic_kind_set)
2780 257714 : atomic_kind => atomic_kind_set(j)
2781 257714 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
2782 257714 : found = .FALSE.
2783 257714 : found_a = .FALSE.
2784 257714 : found_b = .FALSE.
2785 257714 : name_atm_a = name_atm_a_local
2786 257714 : name_atm_b = name_atm_b_local
2787 257714 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
2788 257714 : CALL uppercase(name_atm_a)
2789 257714 : CALL uppercase(name_atm_b)
2790 257714 : pot => potparm_nonbond14%pot(i, j)%pot
2791 :
2792 : ! loop over params from GROMOS
2793 257714 : IF (ASSOCIATED(gro_info%nonbond_a_14)) THEN
2794 540 : ii = 0
2795 540 : jj = 0
2796 1800 : DO k = 1, SIZE(gro_info%nonbond_a_14)
2797 1800 : IF (TRIM(name_atm_a) == TRIM(gro_info%nonbond_a_14(k))) THEN
2798 : ii = k
2799 : found_a = .TRUE.
2800 : EXIT
2801 : END IF
2802 : END DO
2803 2364 : DO k = 1, SIZE(gro_info%nonbond_a_14)
2804 2364 : IF (TRIM(name_atm_b) == TRIM(gro_info%nonbond_a_14(k))) THEN
2805 : jj = k
2806 : found_b = .TRUE.
2807 : EXIT
2808 : END IF
2809 : END DO
2810 540 : IF (ii /= 0 .AND. jj /= 0) THEN
2811 540 : CALL pair_potential_lj_create(pot%set(1)%lj)
2812 1080 : pot%type = lj_type
2813 540 : pot%at1 = name_atm_a
2814 540 : pot%at2 = name_atm_b
2815 540 : pot%set(1)%lj%epsilon = 1.0_dp
2816 540 : pot%set(1)%lj%sigma6 = gro_info%nonbond_c6_14(ii, jj)
2817 540 : pot%set(1)%lj%sigma12 = gro_info%nonbond_c12_14(ii, jj)
2818 540 : pot%rcutsq = (10.0_dp*bohr)**2
2819 540 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
2820 540 : found = .TRUE.
2821 : END IF
2822 : END IF
2823 :
2824 : ! Loop over params from CHARMM
2825 257714 : ii = 0
2826 257714 : jj = 0
2827 257714 : IF (ASSOCIATED(chm_info%nonbond_a_14)) THEN
2828 460416 : DO k = 1, SIZE(chm_info%nonbond_a_14)
2829 460416 : IF ((name_atm_a) == (chm_info%nonbond_a_14(k))) THEN
2830 11206 : ii = k
2831 11206 : rmin2_a = chm_info%nonbond_rmin2_14(k)
2832 11206 : epsilon_a = chm_info%nonbond_eps_14(k)
2833 11206 : found_a = .TRUE.
2834 : END IF
2835 : END DO
2836 460416 : DO k = 1, SIZE(chm_info%nonbond_a_14)
2837 460416 : IF ((name_atm_b) == (chm_info%nonbond_a_14(k))) THEN
2838 8888 : jj = k
2839 8888 : rmin2_b = chm_info%nonbond_rmin2_14(k)
2840 8888 : epsilon_b = chm_info%nonbond_eps_14(k)
2841 8888 : found_b = .TRUE.
2842 : END IF
2843 : END DO
2844 : END IF
2845 257714 : IF (ASSOCIATED(chm_info%nonbond_a)) THEN
2846 48317 : IF (.NOT. found_a) THEN
2847 1442209 : DO k = 1, SIZE(chm_info%nonbond_a)
2848 1442209 : IF ((name_atm_a) == (chm_info%nonbond_a(k))) THEN
2849 37045 : ii = k
2850 37045 : rmin2_a = chm_info%nonbond_rmin2(k)
2851 37045 : epsilon_a = chm_info%nonbond_eps(k)
2852 : END IF
2853 : END DO
2854 : END IF
2855 48317 : IF (.NOT. found_b) THEN
2856 1655119 : DO k = 1, SIZE(chm_info%nonbond_a)
2857 1655119 : IF ((name_atm_b) == (chm_info%nonbond_a(k))) THEN
2858 39411 : jj = k
2859 39411 : rmin2_b = chm_info%nonbond_rmin2(k)
2860 39411 : epsilon_b = chm_info%nonbond_eps(k)
2861 : END IF
2862 : END DO
2863 : END IF
2864 : END IF
2865 257714 : IF (ii /= 0 .AND. jj /= 0) THEN
2866 48251 : rmin = rmin2_a + rmin2_b
2867 : ! ABS to allow for mixing the two different sign conventions for epsilon
2868 48251 : epsilon0 = SQRT(ABS(epsilon_a*epsilon_b))
2869 48251 : CALL pair_potential_lj_create(pot%set(1)%lj)
2870 96502 : pot%type = lj_charmm_type
2871 48251 : pot%at1 = name_atm_a
2872 48251 : pot%at2 = name_atm_b
2873 48251 : pot%set(1)%lj%epsilon = epsilon0
2874 48251 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
2875 48251 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
2876 48251 : pot%rcutsq = (10.0_dp*bohr)**2
2877 48251 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
2878 48251 : found = .TRUE.
2879 : END IF
2880 :
2881 : ! Loop over params from AMBER
2882 257714 : IF (ASSOCIATED(amb_info%nonbond_a)) THEN
2883 199334 : ii = 0
2884 199334 : jj = 0
2885 199334 : IF (.NOT. found_a) THEN
2886 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
2887 45258092 : IF ((name_atm_a) == (amb_info%nonbond_a(k))) THEN
2888 199334 : ii = k
2889 199334 : rmin2_a = amb_info%nonbond_rmin2(k)
2890 199334 : epsilon_a = amb_info%nonbond_eps(k)
2891 : END IF
2892 : END DO
2893 : END IF
2894 199334 : IF (.NOT. found_b) THEN
2895 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
2896 45258092 : IF ((name_atm_b) == (amb_info%nonbond_a(k))) THEN
2897 199334 : jj = k
2898 199334 : rmin2_b = amb_info%nonbond_rmin2(k)
2899 199334 : epsilon_b = amb_info%nonbond_eps(k)
2900 : END IF
2901 : END DO
2902 : END IF
2903 199334 : IF (ii /= 0 .AND. jj /= 0) THEN
2904 199334 : rmin = rmin2_a + rmin2_b
2905 : ! ABS to allow for mixing the two different sign conventions for epsilon
2906 199334 : epsilon0 = SQRT(ABS(epsilon_a*epsilon_b))
2907 199334 : CALL pair_potential_lj_create(pot%set(1)%lj)
2908 398668 : pot%type = lj_charmm_type
2909 199334 : pot%at1 = name_atm_a
2910 199334 : pot%at2 = name_atm_b
2911 199334 : pot%set(1)%lj%epsilon = epsilon0
2912 199334 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
2913 199334 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
2914 199334 : pot%rcutsq = (10.0_dp*bohr)**2
2915 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, &
2916 199334 : name_atm_b)
2917 199334 : found = .TRUE.
2918 : END IF
2919 : END IF
2920 :
2921 : ! Always have the input param last to overwrite all the other ones
2922 257714 : IF (ASSOCIATED(inp_info%nonbonded14)) THEN
2923 12120 : DO k = 1, SIZE(inp_info%nonbonded14%pot)
2924 15263 : IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), &
2925 4817 : " with ", TRIM(inp_info%nonbonded14%pot(k)%pot%at1), &
2926 9634 : TRIM(inp_info%nonbonded14%pot(k)%pot%at2)
2927 : IF ((((name_atm_a) == (inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2928 10446 : ((name_atm_b) == (inp_info%nonbonded14%pot(k)%pot%at2))) .OR. &
2929 : (((name_atm_b) == (inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2930 1674 : ((name_atm_a) == (inp_info%nonbonded14%pot(k)%pot%at2)))) THEN
2931 1666 : IF (ff_type%multiple_potential) THEN
2932 0 : CALL pair_potential_single_add(inp_info%nonbonded14%pot(k)%pot, pot)
2933 0 : IF (found) &
2934 : CALL cp_warn(__LOCATION__, &
2935 : "Multiple ONFO declaration: "//TRIM(name_atm_a)// &
2936 0 : " and "//TRIM(name_atm_b)//" ADDING! ")
2937 0 : potparm_nonbond14%pot(i, j)%pot => pot
2938 0 : potparm_nonbond14%pot(j, i)%pot => pot
2939 : ELSE
2940 1666 : CALL pair_potential_single_copy(inp_info%nonbonded14%pot(k)%pot, pot)
2941 1666 : IF (found) &
2942 : CALL cp_warn(__LOCATION__, &
2943 : "Multiple ONFO declarations: "//TRIM(name_atm_a)// &
2944 0 : " and "//TRIM(name_atm_b)//" OVERWRITING! ")
2945 : END IF
2946 1666 : IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
2947 1666 : found = .TRUE.
2948 : END IF
2949 : END DO
2950 : END IF
2951 :
2952 : ! At the very end we offer the possibility to overwrite the parameters for QM/MM
2953 : ! nonbonded interactions
2954 257714 : IF (use_qmmm_ff) THEN
2955 252 : match_names = 0
2956 252 : IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1
2957 252 : IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1
2958 252 : IF (match_names == 1) THEN
2959 102 : IF (ASSOCIATED(qmmm_env%inp_info%nonbonded14)) THEN
2960 0 : DO k = 1, SIZE(qmmm_env%inp_info%nonbonded14%pot)
2961 : IF (debug_this_module) THEN
2962 : IF (iw > 0) THEN
2963 : WRITE (UNIT=iw, FMT="(T2,A)") &
2964 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
2965 : " with "//TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)//"-"// &
2966 : TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2)
2967 : END IF
2968 : END IF
2969 : IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2970 0 : ((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2))) .OR. &
2971 : (((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2972 0 : ((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2)))) THEN
2973 0 : IF (qmmm_env%multiple_potential) THEN
2974 0 : CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded14%pot(k)%pot, pot)
2975 0 : IF (found) &
2976 : CALL cp_warn(__LOCATION__, &
2977 : "Multiple ONFO declaration: "//TRIM(name_atm_a)// &
2978 0 : " and "//TRIM(name_atm_b)//" Adding QM/MM forcefield specifications")
2979 0 : potparm_nonbond14%pot(i, j)%pot => pot
2980 0 : potparm_nonbond14%pot(j, i)%pot => pot
2981 : ELSE
2982 0 : CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded14%pot(k)%pot, pot)
2983 0 : IF (found) &
2984 : CALL cp_warn(__LOCATION__, &
2985 : "Multiple ONFO declaration: "//TRIM(name_atm_a)// &
2986 0 : " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! ")
2987 : END IF
2988 0 : IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), &
2989 0 : " ", TRIM(name_atm_b)
2990 0 : found = .TRUE.
2991 : END IF
2992 : END DO
2993 : END IF
2994 : END IF
2995 : END IF
2996 :
2997 257714 : IF (.NOT. found) THEN
2998 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
2999 : atm2=TRIM(name_atm_b), &
3000 : type_name="Spline_Bond_Env", &
3001 7923 : array=Ainfo)
3002 7923 : CALL pair_potential_single_clean(pot)
3003 15846 : pot%type = nn_type
3004 7923 : pot%at1 = name_atm_a
3005 7923 : pot%at2 = name_atm_b
3006 : END IF
3007 :
3008 : ! If defined global RCUT let's use it
3009 257714 : IF (ff_type%rcut_nb > 0.0_dp) THEN
3010 26946 : pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb
3011 : END IF
3012 :
3013 : ! Cutoff is defined always as the maximum between the FF and Ewald
3014 257714 : pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut)
3015 268919 : IF (only_qm) THEN
3016 11786 : CALL pair_potential_single_clean(pot)
3017 : END IF
3018 :
3019 : END DO ! atom kind j
3020 :
3021 : END DO ! atom kind i
3022 :
3023 2627 : CALL timestop(handle2)
3024 :
3025 2627 : END SUBROUTINE force_field_pack_nonbond14
3026 :
3027 : ! **************************************************************************************************
3028 : !> \brief Assign input and potential info to potparm_nonbond
3029 : !> \param atomic_kind_set ...
3030 : !> \param ff_type ...
3031 : !> \param qmmm_env ...
3032 : !> \param fatal ...
3033 : !> \param iw ...
3034 : !> \param Ainfo ...
3035 : !> \param chm_info ...
3036 : !> \param inp_info ...
3037 : !> \param gro_info ...
3038 : !> \param amb_info ...
3039 : !> \param potparm_nonbond ...
3040 : !> \param ewald_env ...
3041 : ! **************************************************************************************************
3042 2627 : SUBROUTINE force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, fatal, &
3043 : iw, Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond, &
3044 : ewald_env)
3045 :
3046 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3047 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
3048 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
3049 : LOGICAL :: fatal
3050 : INTEGER :: iw
3051 : CHARACTER(LEN=default_string_length), &
3052 : DIMENSION(:), POINTER :: Ainfo
3053 : TYPE(charmm_info_type), POINTER :: chm_info
3054 : TYPE(input_info_type), POINTER :: inp_info
3055 : TYPE(gromos_info_type), POINTER :: gro_info
3056 : TYPE(amber_info_type), POINTER :: amb_info
3057 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond
3058 : TYPE(ewald_environment_type), POINTER :: ewald_env
3059 :
3060 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond'
3061 :
3062 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
3063 : name_atm_b, name_atm_b_local
3064 : INTEGER :: handle2, i, ii, j, jj, k, match_names
3065 : LOGICAL :: found, is_a_shell, is_b_shell, only_qm, &
3066 : use_qmmm_ff
3067 : REAL(KIND=dp) :: epsilon0, ewald_rcut, rmin
3068 : TYPE(atomic_kind_type), POINTER :: atomic_kind
3069 : TYPE(pair_potential_single_type), POINTER :: pot
3070 :
3071 2627 : CALL timeset(routineN, handle2)
3072 :
3073 2627 : use_qmmm_ff = qmmm_env%use_qmmm_ff
3074 2627 : NULLIFY (pot)
3075 :
3076 2627 : IF (iw > 0) THEN
3077 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
3078 238 : "FORCEFIELD| Checking for nonbonded terms"
3079 : END IF
3080 :
3081 2627 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
3082 2627 : CALL pair_potential_pp_create(potparm_nonbond, SIZE(atomic_kind_set))
3083 :
3084 13832 : DO i = 1, SIZE(atomic_kind_set)
3085 :
3086 11205 : atomic_kind => atomic_kind_set(i)
3087 :
3088 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local, &
3089 11205 : shell_active=is_a_shell)
3090 :
3091 271546 : DO j = i, SIZE(atomic_kind_set)
3092 :
3093 257714 : atomic_kind => atomic_kind_set(j)
3094 :
3095 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local, &
3096 257714 : shell_active=is_b_shell)
3097 :
3098 257714 : found = .FALSE.
3099 :
3100 257714 : name_atm_a = name_atm_a_local
3101 257714 : name_atm_b = name_atm_b_local
3102 257714 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
3103 257714 : CALL uppercase(name_atm_a)
3104 257714 : CALL uppercase(name_atm_b)
3105 257714 : pot => potparm_nonbond%pot(i, j)%pot
3106 :
3107 257714 : IF (iw > 0) THEN
3108 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
3109 : "FORCEFIELD| Checking for nonbonded terms between the atomic kinds "// &
3110 6441 : TRIM(name_atm_a)//" and "//TRIM(name_atm_b)
3111 : END IF
3112 :
3113 : ! Loop over params from GROMOS
3114 257714 : IF (ASSOCIATED(gro_info%nonbond_a)) THEN
3115 540 : ii = 0
3116 540 : jj = 0
3117 1800 : DO k = 1, SIZE(gro_info%nonbond_a)
3118 1800 : IF (TRIM(name_atm_a) == TRIM(gro_info%nonbond_a(k))) THEN
3119 : ii = k
3120 : EXIT
3121 : END IF
3122 : END DO
3123 2364 : DO k = 1, SIZE(gro_info%nonbond_a)
3124 2364 : IF (TRIM(name_atm_b) == TRIM(gro_info%nonbond_a(k))) THEN
3125 : jj = k
3126 : EXIT
3127 : END IF
3128 : END DO
3129 :
3130 540 : IF (ii /= 0 .AND. jj /= 0) THEN
3131 540 : CALL pair_potential_lj_create(pot%set(1)%lj)
3132 1080 : pot%type = lj_type
3133 540 : pot%at1 = name_atm_a
3134 540 : pot%at2 = name_atm_b
3135 540 : pot%set(1)%lj%epsilon = 1.0_dp
3136 540 : pot%set(1)%lj%sigma6 = gro_info%nonbond_c6(ii, jj)
3137 540 : pot%set(1)%lj%sigma12 = gro_info%nonbond_c12(ii, jj)
3138 540 : pot%rcutsq = (10.0_dp*bohr)**2
3139 540 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
3140 540 : found = .TRUE.
3141 : END IF
3142 : END IF
3143 :
3144 : ! Loop over params from CHARMM
3145 257714 : IF (ASSOCIATED(chm_info%nonbond_a)) THEN
3146 48317 : ii = 0
3147 48317 : jj = 0
3148 2286521 : DO k = 1, SIZE(chm_info%nonbond_a)
3149 2286521 : IF ((name_atm_a) == (chm_info%nonbond_a(k))) THEN
3150 48251 : ii = k
3151 : END IF
3152 : END DO
3153 2286521 : DO k = 1, SIZE(chm_info%nonbond_a)
3154 2286521 : IF ((name_atm_b) == (chm_info%nonbond_a(k))) THEN
3155 48299 : jj = k
3156 : END IF
3157 : END DO
3158 :
3159 48317 : IF (ii /= 0 .AND. jj /= 0) THEN
3160 48251 : rmin = chm_info%nonbond_rmin2(ii) + chm_info%nonbond_rmin2(jj)
3161 : epsilon0 = SQRT(chm_info%nonbond_eps(ii)* &
3162 48251 : chm_info%nonbond_eps(jj))
3163 48251 : CALL pair_potential_lj_create(pot%set(1)%lj)
3164 96502 : pot%type = lj_charmm_type
3165 48251 : pot%at1 = name_atm_a
3166 48251 : pot%at2 = name_atm_b
3167 48251 : pot%set(1)%lj%epsilon = epsilon0
3168 48251 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
3169 48251 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
3170 48251 : pot%rcutsq = (10.0_dp*bohr)**2
3171 48251 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
3172 48251 : found = .TRUE.
3173 : END IF
3174 : END IF
3175 :
3176 : ! Loop over params from AMBER
3177 257714 : IF (ASSOCIATED(amb_info%nonbond_a)) THEN
3178 199334 : ii = 0
3179 199334 : jj = 0
3180 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
3181 45258092 : IF ((name_atm_a) == (amb_info%nonbond_a(k))) THEN
3182 199334 : ii = k
3183 : END IF
3184 : END DO
3185 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
3186 45258092 : IF ((name_atm_b) == (amb_info%nonbond_a(k))) THEN
3187 199334 : jj = k
3188 : END IF
3189 : END DO
3190 :
3191 199334 : IF (ii /= 0 .AND. jj /= 0) THEN
3192 199334 : rmin = amb_info%nonbond_rmin2(ii) + amb_info%nonbond_rmin2(jj)
3193 199334 : epsilon0 = SQRT(amb_info%nonbond_eps(ii)*amb_info%nonbond_eps(jj))
3194 199334 : CALL pair_potential_lj_create(pot%set(1)%lj)
3195 398668 : pot%type = lj_charmm_type
3196 199334 : pot%at1 = name_atm_a
3197 199334 : pot%at2 = name_atm_b
3198 199334 : pot%set(1)%lj%epsilon = epsilon0
3199 199334 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
3200 199334 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
3201 199334 : pot%rcutsq = (10.0_dp*bohr)**2
3202 199334 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
3203 199334 : found = .TRUE.
3204 : END IF
3205 : END IF
3206 :
3207 : ! Always have the input param last to overwrite all the other ones
3208 257714 : IF (ASSOCIATED(inp_info%nonbonded)) THEN
3209 53230 : DO k = 1, SIZE(inp_info%nonbonded%pot)
3210 43561 : IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .OR. &
3211 : (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE
3212 : IF (debug_this_module) THEN
3213 : IF (iw > 0) THEN
3214 : WRITE (UNIT=iw, FMT="(T2,A)") &
3215 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3216 : " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// &
3217 : TRIM(inp_info%nonbonded%pot(k)%pot%at2)
3218 : END IF
3219 : END IF
3220 : IF ((((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3221 43559 : ((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at2))) .OR. &
3222 : (((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3223 9669 : ((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at2)))) THEN
3224 9448 : IF (iw > 0) THEN
3225 : WRITE (UNIT=iw, FMT="(T2,A)") &
3226 : "FORCEFIELD| Found nonbonded term "// &
3227 932 : TRIM(name_atm_a)//"-"//TRIM(name_atm_b)
3228 : END IF
3229 9448 : IF (ff_type%multiple_potential) THEN
3230 38 : CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot)
3231 38 : IF (found) &
3232 : CALL cp_warn(__LOCATION__, &
3233 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3234 8 : "-"//TRIM(name_atm_b)//" -> ADDING")
3235 38 : potparm_nonbond%pot(i, j)%pot => pot
3236 38 : potparm_nonbond%pot(j, i)%pot => pot
3237 : ELSE
3238 9410 : CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot)
3239 9410 : IF (found) &
3240 : CALL cp_warn(__LOCATION__, &
3241 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3242 8 : "-"//TRIM(name_atm_b)//" -> OVERWRITING")
3243 : END IF
3244 9448 : found = .TRUE.
3245 : END IF
3246 : END DO
3247 :
3248 : ! Check for wildcards for one of the two types (if not associated yet)
3249 9669 : IF (.NOT. found) THEN
3250 590 : DO k = 1, SIZE(inp_info%nonbonded%pot)
3251 433 : IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .EQV. &
3252 : (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE
3253 : IF (debug_this_module) THEN
3254 : IF (iw > 0) THEN
3255 : WRITE (UNIT=iw, FMT="(T2,A)") &
3256 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3257 : " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// &
3258 : TRIM(inp_info%nonbonded%pot(k)%pot%at2)
3259 : END IF
3260 : END IF
3261 : IF ((name_atm_a == inp_info%nonbonded%pot(k)%pot%at1) .OR. &
3262 : (name_atm_b == inp_info%nonbonded%pot(k)%pot%at2) .OR. &
3263 0 : (name_atm_b == inp_info%nonbonded%pot(k)%pot%at1) .OR. &
3264 157 : (name_atm_a == inp_info%nonbonded%pot(k)%pot%at2)) THEN
3265 0 : IF (iw > 0) THEN
3266 : WRITE (UNIT=iw, FMT="(T2,A)") &
3267 : "FORCEFIELD| Found one wildcard for "// &
3268 0 : TRIM(name_atm_a)//"-"//TRIM(name_atm_b)
3269 : END IF
3270 0 : IF (ff_type%multiple_potential) THEN
3271 0 : CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot)
3272 0 : IF (found) &
3273 : CALL cp_warn(__LOCATION__, &
3274 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3275 0 : "-"//TRIM(name_atm_b)//" -> ADDING")
3276 0 : potparm_nonbond%pot(i, j)%pot => pot
3277 0 : potparm_nonbond%pot(j, i)%pot => pot
3278 : ELSE
3279 0 : CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot)
3280 0 : IF (found) &
3281 : CALL cp_warn(__LOCATION__, &
3282 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3283 0 : "-"//TRIM(name_atm_b)//" -> OVERWRITING")
3284 : END IF
3285 0 : found = .TRUE.
3286 : END IF
3287 : END DO
3288 : END IF
3289 :
3290 : ! Check for wildcards for both types (if not associated yet)
3291 9669 : IF (.NOT. found) THEN
3292 590 : DO k = 1, SIZE(inp_info%nonbonded%pot)
3293 433 : IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) /= "*") .OR. &
3294 : (TRIM(inp_info%nonbonded%pot(k)%pot%at2) /= "*")) CYCLE
3295 : IF (debug_this_module) THEN
3296 : IF (iw > 0) THEN
3297 : WRITE (UNIT=iw, FMT="(T2,A)") &
3298 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3299 : " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// &
3300 : TRIM(inp_info%nonbonded%pot(k)%pot%at2)
3301 : END IF
3302 : END IF
3303 2 : IF (iw > 0) THEN
3304 : WRITE (UNIT=iw, FMT="(T2,A)") &
3305 : "FORCEFIELD| Found wildcards for both "// &
3306 0 : TRIM(name_atm_a)//" and "//TRIM(name_atm_b)
3307 : END IF
3308 2 : IF (ff_type%multiple_potential) THEN
3309 0 : CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot)
3310 0 : IF (found) &
3311 : CALL cp_warn(__LOCATION__, &
3312 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3313 0 : " - "//TRIM(name_atm_b)//" -> ADDING")
3314 0 : potparm_nonbond%pot(i, j)%pot => pot
3315 0 : potparm_nonbond%pot(j, i)%pot => pot
3316 : ELSE
3317 2 : CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot)
3318 2 : IF (found) &
3319 : CALL cp_warn(__LOCATION__, &
3320 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3321 0 : " - "//TRIM(name_atm_b)//" -> OVERWRITING")
3322 : END IF
3323 590 : found = .TRUE.
3324 : END DO
3325 : END IF
3326 : END IF
3327 :
3328 : ! At the very end we offer the possibility to overwrite the parameters for QM/MM
3329 : ! nonbonded interactions
3330 257714 : IF (use_qmmm_ff) THEN
3331 252 : match_names = 0
3332 252 : IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1
3333 252 : IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1
3334 252 : IF (match_names == 1) THEN
3335 102 : IF (ASSOCIATED(qmmm_env%inp_info%nonbonded)) THEN
3336 276 : DO k = 1, SIZE(qmmm_env%inp_info%nonbonded%pot)
3337 : IF (debug_this_module) THEN
3338 : IF (iw > 0) THEN
3339 : WRITE (UNIT=iw, FMT="(T2,A)") &
3340 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3341 : " with "//TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1), &
3342 : TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2)
3343 : END IF
3344 : END IF
3345 : IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3346 174 : ((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2))) .OR. &
3347 : (((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3348 102 : ((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2)))) THEN
3349 20 : IF (iw > 0) THEN
3350 : WRITE (UNIT=iw, FMT="(T2,A)") &
3351 11 : "FORCEFIELD| Found "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)//" (QM/MM)"
3352 : END IF
3353 20 : IF (qmmm_env%multiple_potential) THEN
3354 0 : CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot)
3355 0 : IF (found) &
3356 : CALL cp_warn(__LOCATION__, &
3357 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3358 0 : " and "//TRIM(name_atm_b)//" -> ADDING QM/MM forcefield specifications")
3359 0 : potparm_nonbond%pot(i, j)%pot => pot
3360 0 : potparm_nonbond%pot(j, i)%pot => pot
3361 : ELSE
3362 20 : CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot)
3363 20 : IF (found) &
3364 : CALL cp_warn(__LOCATION__, &
3365 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3366 2 : " and "//TRIM(name_atm_b)//" -> OVERWRITING QM/MM forcefield specifications")
3367 : END IF
3368 20 : found = .TRUE.
3369 : END IF
3370 : END DO
3371 : END IF
3372 : END IF
3373 : END IF
3374 :
3375 257714 : IF (.NOT. found) THEN
3376 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
3377 : atm2=TRIM(name_atm_b), &
3378 : type_name="Spline_Non_Bond_Env", &
3379 : fatal=fatal, &
3380 137 : array=Ainfo)
3381 : END IF
3382 :
3383 : ! If defined global RCUT let's use it
3384 257714 : IF (ff_type%rcut_nb > 0.0_dp) THEN
3385 26946 : pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb
3386 : END IF
3387 :
3388 : ! Cutoff is defined always as the maximum between the FF and Ewald
3389 257714 : pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut)
3390 : ! Set the shell type
3391 257714 : IF ((is_a_shell .AND. .NOT. is_b_shell) .OR. (is_b_shell .AND. .NOT. is_a_shell)) THEN
3392 66 : pot%shell_type = nosh_sh
3393 257648 : ELSE IF (is_a_shell .AND. is_b_shell) THEN
3394 616 : pot%shell_type = sh_sh
3395 : ELSE
3396 257032 : pot%shell_type = nosh_nosh
3397 : END IF
3398 :
3399 526633 : IF (only_qm) THEN
3400 11786 : CALL pair_potential_single_clean(pot)
3401 : END IF
3402 :
3403 : END DO ! jkind
3404 :
3405 : END DO ! ikind
3406 :
3407 2627 : CALL timestop(handle2)
3408 :
3409 2627 : END SUBROUTINE force_field_pack_nonbond
3410 :
3411 : ! **************************************************************************************************
3412 : !> \brief create the pair potential spline environment
3413 : !> \param atomic_kind_set ...
3414 : !> \param ff_type ...
3415 : !> \param iw2 ...
3416 : !> \param iw3 ...
3417 : !> \param iw4 ...
3418 : !> \param potparm ...
3419 : !> \param do_zbl ...
3420 : !> \param nonbonded_type ...
3421 : ! **************************************************************************************************
3422 5254 : SUBROUTINE force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, &
3423 : potparm, do_zbl, nonbonded_type)
3424 :
3425 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3426 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
3427 : INTEGER :: iw2, iw3, iw4
3428 : TYPE(pair_potential_pp_type), POINTER :: potparm
3429 : LOGICAL, INTENT(IN) :: do_zbl
3430 : CHARACTER(LEN=*), INTENT(IN) :: nonbonded_type
3431 :
3432 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_splines'
3433 :
3434 : INTEGER :: handle2, ikind, jkind, n
3435 5254 : TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spl_p
3436 : TYPE(spline_environment_type), POINTER :: spline_env
3437 :
3438 5254 : CALL timeset(routineN, handle2)
3439 :
3440 5254 : IF (iw2 > 0) THEN
3441 : WRITE (UNIT=iw2, FMT="(/,T2,A)") &
3442 476 : "FORCEFIELD| Splining nonbonded terms"
3443 : END IF
3444 :
3445 : ! Figure out which nonbonded interactions happen to be identical, and
3446 : ! prepare storage for these, avoiding duplicates.
3447 5254 : NULLIFY (spline_env)
3448 : CALL get_nonbond_storage(spline_env, potparm, atomic_kind_set, &
3449 5254 : do_zbl, shift_cutoff=ff_type%shift_cutoff)
3450 : ! Effectively compute the spline data
3451 : CALL spline_nonbond_control(spline_env, potparm, &
3452 : atomic_kind_set, eps_spline=ff_type%eps_spline, &
3453 : max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb, &
3454 : emax_spline=ff_type%emax_spline, npoints=ff_type%npoints, &
3455 : iw=iw2, iw2=iw3, iw3=iw4, &
3456 : do_zbl=do_zbl, shift_cutoff=ff_type%shift_cutoff, &
3457 5254 : nonbonded_type=nonbonded_type)
3458 : ! Let the pointers on potparm point to the splines generated in
3459 : ! spline_nonbond_control
3460 27664 : DO ikind = 1, SIZE(potparm%pot, 1)
3461 543092 : DO jkind = ikind, SIZE(potparm%pot, 2)
3462 515428 : n = spline_env%spltab(ikind, jkind)
3463 515428 : spl_p => spline_env%spl_pp(n)%spl_p
3464 515428 : CALL spline_data_p_retain(spl_p)
3465 515428 : CALL spline_data_p_release(potparm%pot(ikind, jkind)%pot%pair_spline_data)
3466 537838 : potparm%pot(ikind, jkind)%pot%pair_spline_data => spl_p
3467 : END DO
3468 : END DO
3469 5254 : CALL spline_env_release(spline_env)
3470 5254 : DEALLOCATE (spline_env)
3471 : NULLIFY (spline_env)
3472 :
3473 5254 : IF (iw2 > 0) THEN
3474 : WRITE (UNIT=iw2, FMT="(/,T2,A)") &
3475 476 : "FORCEFIELD| Splining done"
3476 : END IF
3477 :
3478 5254 : CALL timestop(handle2)
3479 :
3480 5254 : END SUBROUTINE force_field_pack_splines
3481 :
3482 : ! **************************************************************************************************
3483 : !> \brief Compute the electrostatic interaction cutoffs
3484 : !> \param atomic_kind_set ...
3485 : !> \param ff_type ...
3486 : !> \param potparm_nonbond ...
3487 : !> \param ewald_env ...
3488 : !> \param iw ...
3489 : !> \author Toon.Verstraelen@gmail.com
3490 : ! **************************************************************************************************
3491 2643 : SUBROUTINE force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, ewald_env, iw)
3492 :
3493 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3494 : TYPE(force_field_type), INTENT(IN) :: ff_type
3495 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond
3496 : TYPE(ewald_environment_type), POINTER :: ewald_env
3497 : INTEGER, INTENT(IN) :: iw
3498 :
3499 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_eicut'
3500 :
3501 : INTEGER :: ewald_type, handle, i1, i2, nkinds
3502 : REAL(KIND=dp) :: alpha, beta, mm_radius1, mm_radius2, &
3503 : rcut2, rcut2_ewald, tmp
3504 2643 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs
3505 : TYPE(atomic_kind_type), POINTER :: atomic_kind
3506 :
3507 2643 : CALL timeset(routineN, handle)
3508 :
3509 2643 : IF (iw > 0) THEN
3510 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
3511 239 : "FORCEFIELD| Computing the electrostatic interactions cutoffs"
3512 : END IF
3513 :
3514 2643 : tmp = 0.0_dp
3515 2643 : nkinds = SIZE(atomic_kind_set)
3516 :
3517 : ! Allocate the array with interaction cutoffs for the electrostatics, used
3518 : ! to make the electrostatic interaction continuous at ewald_env%rcut
3519 10572 : ALLOCATE (interaction_cutoffs(3, nkinds, nkinds))
3520 2032934 : interaction_cutoffs = 0.0_dp
3521 :
3522 : ! Compute the interaction cutoff if SHIFT_CUTOFF is active
3523 2643 : IF (ff_type%shift_cutoff) THEN
3524 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type, &
3525 2491 : rcut=rcut2_ewald)
3526 2491 : rcut2_ewald = rcut2_ewald*rcut2_ewald
3527 11730 : DO i1 = 1, nkinds
3528 9239 : atomic_kind => atomic_kind_set(i1)
3529 9239 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius1)
3530 118103 : DO i2 = 1, nkinds
3531 106373 : rcut2 = rcut2_ewald
3532 106373 : IF (ASSOCIATED(potparm_nonbond)) THEN
3533 105843 : rcut2 = MAX(potparm_nonbond%pot(i1, i2)%pot%rcutsq, rcut2_ewald)
3534 : END IF
3535 115612 : IF (rcut2 > 0) THEN
3536 103061 : atomic_kind => atomic_kind_set(i2)
3537 103061 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius2)
3538 : ! cutoff for core-core
3539 : interaction_cutoffs(1, i1, i2) = potential_coulomb(rcut2, tmp, &
3540 103061 : 1.0_dp, ewald_type, alpha, 0.0_dp, 0.0_dp)
3541 : ! cutoff for core-shell, core-ion, shell-core or ion-core
3542 103061 : IF (mm_radius1 > 0.0_dp) THEN
3543 676 : beta = sqrthalf/mm_radius1
3544 : ELSE
3545 102385 : beta = 0.0_dp
3546 : END IF
3547 : interaction_cutoffs(2, i1, i2) = potential_coulomb(rcut2, tmp, &
3548 103061 : 1.0_dp, ewald_type, alpha, beta, 0.0_dp)
3549 : ! cutoff for shell-shell or ion-ion
3550 103061 : IF (mm_radius1 + mm_radius2 > 0.0_dp) THEN
3551 698 : beta = sqrthalf/SQRT(mm_radius1*mm_radius1 + mm_radius2*mm_radius2)
3552 : ELSE
3553 102363 : beta = 0.0_dp
3554 : END IF
3555 : interaction_cutoffs(3, i1, i2) = potential_coulomb(rcut2, tmp, &
3556 103061 : 1.0_dp, ewald_type, alpha, beta, 0.0_dp)
3557 : END IF
3558 : END DO
3559 : END DO
3560 : END IF
3561 :
3562 2643 : CALL ewald_env_set(ewald_env, interaction_cutoffs=interaction_cutoffs)
3563 :
3564 2643 : CALL timestop(handle)
3565 :
3566 2643 : END SUBROUTINE force_field_pack_eicut
3567 :
3568 : ! **************************************************************************************************
3569 : !> \brief Issues on screen a warning when repetitions are present in the
3570 : !> definition of the forcefield
3571 : !> \param found ...
3572 : !> \param tag_label ...
3573 : !> \param name_atm_a ...
3574 : !> \param name_atm_b ...
3575 : !> \param name_atm_c ...
3576 : !> \param name_atm_d ...
3577 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
3578 : ! **************************************************************************************************
3579 781207 : SUBROUTINE issue_duplications(found, tag_label, name_atm_a, name_atm_b, &
3580 : name_atm_c, name_atm_d)
3581 :
3582 : LOGICAL, INTENT(IN) :: found
3583 : CHARACTER(LEN=*), INTENT(IN) :: tag_label, name_atm_a
3584 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name_atm_b, name_atm_c, name_atm_d
3585 :
3586 : CHARACTER(LEN=default_string_length) :: item
3587 :
3588 781207 : item = "("//TRIM(name_atm_a)
3589 781207 : IF (PRESENT(name_atm_b)) THEN
3590 773361 : item = TRIM(item)//", "//TRIM(name_atm_b)
3591 : END IF
3592 781207 : IF (PRESENT(name_atm_c)) THEN
3593 163056 : item = TRIM(item)//", "//TRIM(name_atm_c)
3594 : END IF
3595 781207 : IF (PRESENT(name_atm_d)) THEN
3596 3418 : item = TRIM(item)//", "//TRIM(name_atm_d)
3597 : END IF
3598 781207 : item = TRIM(item)//")"
3599 781207 : IF (found) THEN
3600 1678 : CPWARN("Found multiple "//TRIM(tag_label)//" terms for "//TRIM(item)//" -> OVERWRITING")
3601 : END IF
3602 :
3603 781207 : END SUBROUTINE issue_duplications
3604 :
3605 : ! **************************************************************************************************
3606 : !> \brief Store informations on possible missing ForceFields parameters
3607 : !> \param atm1 ...
3608 : !> \param atm2 ...
3609 : !> \param atm3 ...
3610 : !> \param atm4 ...
3611 : !> \param type_name ...
3612 : !> \param fatal ...
3613 : !> \param array ...
3614 : ! **************************************************************************************************
3615 167316 : SUBROUTINE store_FF_missing_par(atm1, atm2, atm3, atm4, type_name, fatal, array)
3616 : CHARACTER(LEN=*), INTENT(IN) :: atm1
3617 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: atm2, atm3, atm4
3618 : CHARACTER(LEN=*), INTENT(IN) :: type_name
3619 : LOGICAL, INTENT(INOUT), OPTIONAL :: fatal
3620 : CHARACTER(LEN=default_string_length), &
3621 : DIMENSION(:), POINTER :: array
3622 :
3623 : CHARACTER(LEN=10) :: sfmt
3624 : CHARACTER(LEN=9) :: my_atm1, my_atm2, my_atm3, my_atm4
3625 : CHARACTER(LEN=default_path_length) :: my_format
3626 : INTEGER :: fmt, i, nsize
3627 : LOGICAL :: found
3628 :
3629 167316 : nsize = 0
3630 167316 : fmt = 1
3631 : my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)// &
3632 167316 : '",T40,"(",A9,")")'
3633 167316 : IF (PRESENT(atm2)) fmt = fmt + 1
3634 167316 : IF (PRESENT(atm3)) fmt = fmt + 1
3635 167316 : IF (PRESENT(atm4)) fmt = fmt + 1
3636 167316 : CALL integer_to_string(fmt - 1, sfmt)
3637 167316 : IF (fmt > 1) &
3638 : my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)// &
3639 167304 : '",T40,"(",A9,'//TRIM(sfmt)//'(",",A9),")")'
3640 167316 : IF (PRESENT(fatal)) fatal = .TRUE.
3641 : ! Check for previous already stored equal force fields
3642 167316 : IF (ASSOCIATED(array)) nsize = SIZE(array)
3643 167316 : found = .FALSE.
3644 167316 : IF (nsize >= 1) THEN
3645 19478780 : DO i = 1, nsize
3646 8 : SELECT CASE (type_name)
3647 : CASE ("Bond")
3648 8 : IF (INDEX(array(i) (21:39), "Bond") == 0) CYCLE
3649 8 : my_atm1 = array(i) (41:49)
3650 8 : my_atm2 = array(i) (51:59)
3651 8 : CALL compress(my_atm1, .TRUE.)
3652 8 : CALL compress(my_atm2, .TRUE.)
3653 8 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. &
3654 8 : ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .TRUE.
3655 : CASE ("Angle")
3656 8 : IF (INDEX(array(i) (21:39), "Angle") == 0) CYCLE
3657 0 : my_atm1 = array(i) (41:49)
3658 0 : my_atm2 = array(i) (51:59)
3659 0 : my_atm3 = array(i) (61:69)
3660 0 : CALL compress(my_atm1, .TRUE.)
3661 0 : CALL compress(my_atm2, .TRUE.)
3662 0 : CALL compress(my_atm3, .TRUE.)
3663 0 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. &
3664 : ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) &
3665 18203730 : found = .TRUE.
3666 : CASE ("Urey-Bradley")
3667 18203730 : IF (INDEX(array(i) (21:39), "Urey-Bradley") == 0) CYCLE
3668 18203730 : my_atm1 = array(i) (41:49)
3669 18203730 : my_atm2 = array(i) (51:59)
3670 18203730 : my_atm3 = array(i) (61:69)
3671 18203730 : CALL compress(my_atm1, .TRUE.)
3672 18203730 : CALL compress(my_atm2, .TRUE.)
3673 18203730 : CALL compress(my_atm3, .TRUE.)
3674 18203730 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. &
3675 : ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) &
3676 600240 : found = .TRUE.
3677 : CASE ("Torsion")
3678 600240 : IF (INDEX(array(i) (21:39), "Torsion") == 0) CYCLE
3679 196462 : my_atm1 = array(i) (41:49)
3680 196462 : my_atm2 = array(i) (51:59)
3681 196462 : my_atm3 = array(i) (61:69)
3682 196462 : my_atm4 = array(i) (71:79)
3683 196462 : CALL compress(my_atm1, .TRUE.)
3684 196462 : CALL compress(my_atm2, .TRUE.)
3685 196462 : CALL compress(my_atm3, .TRUE.)
3686 196462 : CALL compress(my_atm4, .TRUE.)
3687 196462 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3688 : ((atm1 == my_atm4) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm1))) &
3689 154212 : found = .TRUE.
3690 : CASE ("Improper")
3691 154212 : IF (INDEX(array(i) (21:39), "Improper") == 0) CYCLE
3692 9684 : my_atm1 = array(i) (41:49)
3693 9684 : my_atm2 = array(i) (51:59)
3694 9684 : my_atm3 = array(i) (61:69)
3695 9684 : my_atm4 = array(i) (71:79)
3696 9684 : CALL compress(my_atm1, .TRUE.)
3697 9684 : CALL compress(my_atm2, .TRUE.)
3698 9684 : CALL compress(my_atm3, .TRUE.)
3699 9684 : CALL compress(my_atm4, .TRUE.)
3700 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3701 : ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4)) .OR. &
3702 : ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3)) .OR. &
3703 : ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm2)) .OR. &
3704 9684 : ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm3)) .OR. &
3705 : ((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3))) &
3706 483920 : found = .TRUE.
3707 :
3708 : CASE ("Out of plane bend")
3709 483920 : IF (INDEX(array(i) (21:39), "Out of plane bend") == 0) CYCLE
3710 27416 : my_atm1 = array(i) (41:49)
3711 27416 : my_atm2 = array(i) (51:59)
3712 27416 : my_atm3 = array(i) (61:69)
3713 27416 : my_atm4 = array(i) (71:79)
3714 27416 : CALL compress(my_atm1, .TRUE.)
3715 27416 : CALL compress(my_atm2, .TRUE.)
3716 27416 : CALL compress(my_atm3, .TRUE.)
3717 27416 : CALL compress(my_atm4, .TRUE.)
3718 27416 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3719 : ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4))) &
3720 8 : found = .TRUE.
3721 :
3722 : CASE ("Charge")
3723 8 : IF (INDEX(array(i) (21:39), "Charge") == 0) CYCLE
3724 8 : my_atm1 = array(i) (41:49)
3725 8 : CALL compress(my_atm1, .TRUE.)
3726 8 : IF (atm1 == my_atm1) found = .TRUE.
3727 : CASE ("Spline_Bond_Env", "Spline_Non_Bond_Env")
3728 18820 : IF (INDEX(array(i) (21:39), "Spline_") == 0) CYCLE
3729 6567 : fmt = 0
3730 6567 : my_atm1 = array(i) (41:49)
3731 6567 : my_atm2 = array(i) (51:59)
3732 6567 : CALL compress(my_atm1, .TRUE.)
3733 6567 : CALL compress(my_atm2, .TRUE.)
3734 6567 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. &
3735 0 : ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .TRUE.
3736 : CASE DEFAULT
3737 : ! Should never reach this point
3738 19460946 : CPABORT("")
3739 : END SELECT
3740 18315455 : IF (found) EXIT
3741 : END DO
3742 : END IF
3743 167316 : IF (.NOT. found) THEN
3744 21062 : nsize = nsize + 1
3745 21062 : CALL reallocate(array, 1, nsize)
3746 12 : SELECT CASE (fmt)
3747 : CASE (1)
3748 12 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1
3749 : CASE (2)
3750 1501 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2
3751 : CASE (3)
3752 11668 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2, atm3
3753 : CASE (4)
3754 21062 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2, atm3, atm4
3755 : END SELECT
3756 : END IF
3757 :
3758 167316 : END SUBROUTINE store_FF_missing_par
3759 :
3760 : ! **************************************************************************************************
3761 : !> \brief Search sorted 2d array of integers for a first occurence of value `val` in row `row`
3762 : !> \param array 2d array of integers
3763 : !> \param val value to search
3764 : !> \param row row to search, default = 1
3765 : !> \return column index if `val` is found in the row `row` of `array`; zero otherwise
3766 : ! **************************************************************************************************
3767 45098 : FUNCTION bsearch_leftmost_2d(array, val, row) RESULT(res)
3768 : INTEGER, INTENT(IN) :: array(:, :), val
3769 : INTEGER, INTENT(IN), OPTIONAL :: row
3770 : INTEGER :: res
3771 :
3772 : INTEGER :: left, locRow, mid, right
3773 :
3774 45098 : locRow = 1
3775 45098 : IF (PRESENT(row)) locRow = row
3776 :
3777 45098 : left = 1
3778 90196 : right = UBOUND(array, dim=2)
3779 :
3780 571050 : DO WHILE (left < right)
3781 525952 : mid = (left + right)/2
3782 571050 : IF (array(locRow, mid) < val) THEN
3783 349610 : left = mid + 1
3784 : ELSE
3785 : right = mid
3786 : END IF
3787 : END DO
3788 :
3789 45098 : res = left
3790 :
3791 : ! Not found:
3792 45098 : IF (array(locRow, res) /= val) res = 0
3793 :
3794 45098 : END FUNCTION bsearch_leftmost_2d
3795 :
3796 : END MODULE force_fields_all
|
