Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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, mace_type, &
68 : nequip_type, nn_type, nosh_nosh, nosh_sh, pair_potential_lj_create, &
69 : pair_potential_pp_create, pair_potential_pp_type, pair_potential_single_add, &
70 : pair_potential_single_clean, pair_potential_single_copy, pair_potential_single_type, &
71 : sh_sh, siepmann_type, 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
143 2637 : INTEGER, POINTER :: map_bond_kind(:)
144 : LOGICAL :: found
145 : TYPE(atomic_kind_type), POINTER :: atomic_kind
146 2637 : TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set
147 2637 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
148 : TYPE(molecule_kind_type), POINTER :: molecule_kind
149 : TYPE(molecule_type), POINTER :: molecule
150 :
151 2637 : CALL timeset(routineN, handle2)
152 :
153 2637 : IF (iw > 0) THEN
154 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
155 242 : "FORCEFIELD| Checking for unique bond terms"
156 : END IF
157 :
158 74503 : DO i = 1, SIZE(molecule_kind_set)
159 71866 : molecule_kind => molecule_kind_set(i)
160 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
161 : molecule_list=molecule_list, &
162 : natom=natom, &
163 71866 : nbond=nbond, bond_list=bond_list)
164 71866 : molecule => molecule_set(molecule_list(1))
165 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
166 146369 : IF (nbond > 0) THEN
167 88305 : ALLOCATE (map_bond_kind(nbond))
168 29435 : counter = 0
169 29435 : 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 144238 : DO j = 1, nbond
176 114823 : atm_a = bond_list(j)%a
177 114823 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
178 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
179 114823 : name=name_atm_a)
180 114823 : atm_b = bond_list(j)%b
181 114823 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
182 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
183 114823 : name=name_atm_b)
184 114823 : found = .FALSE.
185 483453 : DO k = 1, j - 1
186 416668 : atm_a = bond_list(k)%a
187 416668 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
188 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
189 416668 : name=name_atm_a2)
190 416668 : atm_b = bond_list(k)%b
191 416668 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
192 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
193 416668 : name=name_atm_b2)
194 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
195 416668 : ((name_atm_b) == (name_atm_b2))) .OR. &
196 : (((name_atm_a) == (name_atm_b2)) .AND. &
197 66785 : ((name_atm_b) == (name_atm_a2)))) THEN
198 48038 : found = .TRUE.
199 48038 : map_bond_kind(j) = map_bond_kind(k)
200 : EXIT
201 : END IF
202 : END DO
203 29415 : IF (.NOT. found) THEN
204 66785 : counter = counter + 1
205 66785 : map_bond_kind(j) = counter
206 : END IF
207 : END DO
208 : END IF
209 29435 : NULLIFY (bond_kind_set)
210 29435 : CALL allocate_bond_kind_set(bond_kind_set, counter)
211 144386 : DO j = 1, nbond
212 144386 : 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 29435 : bond_kind_set=bond_kind_set, bond_list=bond_list)
216 29435 : DEALLOCATE (map_bond_kind)
217 : END IF
218 : END DO
219 2637 : CALL timestop(handle2)
220 :
221 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
246 2637 : INTEGER, POINTER :: map_bend_kind(:)
247 : LOGICAL :: found
248 : TYPE(atomic_kind_type), POINTER :: atomic_kind
249 2637 : TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set
250 2637 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
251 : TYPE(molecule_kind_type), POINTER :: molecule_kind
252 : TYPE(molecule_type), POINTER :: molecule
253 :
254 2637 : CALL timeset(routineN, handle2)
255 :
256 2637 : IF (iw > 0) THEN
257 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
258 242 : "FORCEFIELD| Checking for unique bend terms"
259 : END IF
260 :
261 74503 : DO i = 1, SIZE(molecule_kind_set)
262 71866 : molecule_kind => molecule_kind_set(i)
263 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
264 : molecule_list=molecule_list, &
265 : natom=natom, &
266 71866 : nbend=nbend, bend_list=bend_list)
267 71866 : molecule => molecule_set(molecule_list(1))
268 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
269 146369 : IF (nbend > 0) THEN
270 87315 : ALLOCATE (map_bend_kind(nbend))
271 29105 : counter = 0
272 29105 : 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 169571 : DO j = 1, nbend
279 140478 : atm_a = bend_list(j)%a
280 140478 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
281 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
282 140478 : name=name_atm_a)
283 140478 : atm_b = bend_list(j)%b
284 140478 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
285 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
286 140478 : name=name_atm_b)
287 140478 : atm_c = bend_list(j)%c
288 140478 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
289 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
290 140478 : name=name_atm_c)
291 140478 : found = .FALSE.
292 2499741 : DO k = 1, j - 1
293 2404937 : atm_a = bend_list(k)%a
294 2404937 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
295 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
296 2404937 : name=name_atm_a2)
297 2404937 : atm_b = bend_list(k)%b
298 2404937 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
299 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
300 2404937 : name=name_atm_b2)
301 2404937 : atm_c = bend_list(k)%c
302 2404937 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
303 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
304 2404937 : name=name_atm_c2)
305 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
306 : ((name_atm_b) == (name_atm_b2)) .AND. &
307 2404937 : ((name_atm_c) == (name_atm_c2))) .OR. &
308 : (((name_atm_a) == (name_atm_c2)) .AND. &
309 : ((name_atm_b) == (name_atm_b2)) .AND. &
310 94804 : ((name_atm_c) == (name_atm_a2)))) THEN
311 45674 : found = .TRUE.
312 45674 : map_bend_kind(j) = map_bend_kind(k)
313 : EXIT
314 : END IF
315 : END DO
316 29093 : IF (.NOT. found) THEN
317 94804 : counter = counter + 1
318 94804 : map_bend_kind(j) = counter
319 : END IF
320 : END DO
321 : END IF
322 29105 : NULLIFY (bend_kind_set)
323 29105 : CALL allocate_bend_kind_set(bend_kind_set, counter)
324 169739 : DO j = 1, nbend
325 169739 : 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 29105 : bend_kind_set=bend_kind_set, bend_list=bend_list)
329 29105 : DEALLOCATE (map_bend_kind)
330 : END IF
331 : END DO
332 :
333 2637 : CALL timestop(handle2)
334 :
335 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
358 2637 : 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 2637 : TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set
364 2637 : TYPE(ub_type), DIMENSION(:), POINTER :: ub_list
365 :
366 2637 : CALL timeset(routineN, handle2)
367 :
368 2637 : IF (iw > 0) THEN
369 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
370 242 : "FORCEFIELD| Checking for unique Urey-Bradley terms"
371 : END IF
372 :
373 74503 : DO i = 1, SIZE(molecule_kind_set)
374 71866 : molecule_kind => molecule_kind_set(i)
375 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
376 : molecule_list=molecule_list, &
377 : natom=natom, &
378 71866 : nub=nub, ub_list=ub_list)
379 71866 : molecule => molecule_set(molecule_list(1))
380 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
381 146369 : IF (nub > 0) THEN
382 87273 : ALLOCATE (map_ub_kind(nub))
383 29091 : counter = 0
384 169567 : DO j = 1, nub
385 140476 : atm_a = ub_list(j)%a
386 140476 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
387 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
388 140476 : name=name_atm_a)
389 140476 : atm_b = ub_list(j)%b
390 140476 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
391 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
392 140476 : name=name_atm_b)
393 140476 : atm_c = ub_list(j)%c
394 140476 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
395 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
396 140476 : name=name_atm_c)
397 140476 : found = .FALSE.
398 2499739 : DO k = 1, j - 1
399 2404937 : atm_a = ub_list(k)%a
400 2404937 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
401 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
402 2404937 : name=name_atm_a2)
403 2404937 : atm_b = ub_list(k)%b
404 2404937 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
405 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
406 2404937 : name=name_atm_b2)
407 2404937 : atm_c = ub_list(k)%c
408 2404937 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
409 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
410 2404937 : name=name_atm_c2)
411 : IF ((((name_atm_a) == (name_atm_a2)) .AND. &
412 : ((name_atm_b) == (name_atm_b2)) .AND. &
413 2404937 : ((name_atm_c) == (name_atm_c2))) .OR. &
414 : (((name_atm_a) == (name_atm_c2)) .AND. &
415 : ((name_atm_b) == (name_atm_b2)) .AND. &
416 94802 : ((name_atm_c) == (name_atm_a2)))) THEN
417 45674 : found = .TRUE.
418 45674 : map_ub_kind(j) = map_ub_kind(k)
419 : EXIT
420 : END IF
421 : END DO
422 29091 : IF (.NOT. found) THEN
423 94802 : counter = counter + 1
424 94802 : map_ub_kind(j) = counter
425 : END IF
426 : END DO
427 29091 : CALL allocate_ub_kind_set(ub_kind_set, counter)
428 169567 : DO j = 1, nub
429 169567 : 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 29091 : ub_kind_set=ub_kind_set, ub_list=ub_list)
433 29091 : DEALLOCATE (map_ub_kind)
434 : END IF
435 : END DO
436 2637 : CALL timestop(handle2)
437 :
438 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
465 2637 : 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 2637 : TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set
471 2637 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
472 :
473 2637 : CALL timeset(routineN, handle2)
474 :
475 2637 : IF (iw > 0) THEN
476 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
477 242 : "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 2637 : chk_reverse = (ff_type%ff_type /= do_ff_amber)
483 :
484 74503 : DO i = 1, SIZE(molecule_kind_set)
485 71866 : molecule_kind => molecule_kind_set(i)
486 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
487 : molecule_list=molecule_list, &
488 : natom=natom, &
489 71866 : ntorsion=ntorsion, torsion_list=torsion_list)
490 71866 : molecule => molecule_set(molecule_list(1))
491 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
492 146369 : IF (ntorsion > 0) THEN
493 16602 : ALLOCATE (map_torsion_kind(ntorsion))
494 5534 : counter = 0
495 5534 : 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 162887 : DO j = 1, ntorsion
502 157361 : atm_a = torsion_list(j)%a
503 157361 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
504 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
505 157361 : name=name_atm_a)
506 157361 : atm_b = torsion_list(j)%b
507 157361 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
508 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
509 157361 : name=name_atm_b)
510 157361 : atm_c = torsion_list(j)%c
511 157361 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
512 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
513 157361 : name=name_atm_c)
514 157361 : atm_d = torsion_list(j)%d
515 157361 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
516 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
517 157361 : name=name_atm_d)
518 157361 : found = .FALSE.
519 2932946 : DO k = 1, j - 1
520 2840585 : atm_a = torsion_list(k)%a
521 2840585 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
522 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
523 2840585 : name=name_atm_a2)
524 2840585 : atm_b = torsion_list(k)%b
525 2840585 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
526 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
527 2840585 : name=name_atm_b2)
528 2840585 : atm_c = torsion_list(k)%c
529 2840585 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
530 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
531 2840585 : name=name_atm_c2)
532 2840585 : atm_d = torsion_list(k)%d
533 2840585 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
534 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
535 2840585 : 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 2840585 : ((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 92361 : ((name_atm_d) == (name_atm_a2)))) THEN
545 65000 : found = .TRUE.
546 65000 : map_torsion_kind(j) = map_torsion_kind(k)
547 : EXIT
548 : END IF
549 : END DO
550 5526 : IF (.NOT. found) THEN
551 92361 : counter = counter + 1
552 92361 : map_torsion_kind(j) = counter
553 : END IF
554 : END DO
555 : END IF
556 5534 : NULLIFY (torsion_kind_set)
557 5534 : CALL allocate_torsion_kind_set(torsion_kind_set, counter)
558 163207 : DO j = 1, ntorsion
559 163207 : 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 5534 : torsion_kind_set=torsion_kind_set, torsion_list=torsion_list)
563 5534 : DEALLOCATE (map_torsion_kind)
564 : END IF
565 : END DO
566 :
567 2637 : CALL timestop(handle2)
568 :
569 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
596 2637 : INTEGER, POINTER :: map_impr_kind(:)
597 : LOGICAL :: found
598 : TYPE(atomic_kind_type), POINTER :: atomic_kind
599 2637 : TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set
600 2637 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list
601 : TYPE(molecule_kind_type), POINTER :: molecule_kind
602 : TYPE(molecule_type), POINTER :: molecule
603 :
604 2637 : CALL timeset(routineN, handle2)
605 :
606 2637 : IF (iw > 0) THEN
607 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
608 242 : "FORCEFIELD| Checking for unique improper terms"
609 : END IF
610 :
611 74503 : DO i = 1, SIZE(molecule_kind_set)
612 71866 : molecule_kind => molecule_kind_set(i)
613 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
614 : molecule_list=molecule_list, &
615 : natom=natom, &
616 71866 : nimpr=nimpr, impr_list=impr_list)
617 71866 : molecule => molecule_set(molecule_list(1))
618 :
619 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
620 :
621 146369 : 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 2637 : CALL timestop(handle2)
695 :
696 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
725 2637 : 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 2637 : TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set
731 2637 : TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list
732 :
733 2637 : CALL timeset(routineN, handle2)
734 :
735 2637 : IF (iw > 0) THEN
736 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
737 242 : "FORCEFIELD| Checking for unique out-of-plane bend terms"
738 : END IF
739 :
740 74503 : DO i = 1, SIZE(molecule_kind_set)
741 71866 : molecule_kind => molecule_kind_set(i)
742 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
743 : molecule_list=molecule_list, &
744 : natom=natom, &
745 71866 : nopbend=nopbend, opbend_list=opbend_list)
746 71866 : molecule => molecule_set(molecule_list(1))
747 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
748 146369 : 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 2637 : CALL timestop(handle2)
822 :
823 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
859 : LOGICAL :: found, only_qm
860 : TYPE(atomic_kind_type), POINTER :: atomic_kind
861 2637 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
862 : TYPE(molecule_kind_type), POINTER :: molecule_kind
863 : TYPE(molecule_type), POINTER :: molecule
864 :
865 2637 : CALL timeset(routineN, handle2)
866 :
867 2637 : IF (iw > 0) THEN
868 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
869 242 : "FORCEFIELD| Checking for bond terms"
870 : END IF
871 :
872 74503 : DO i = 1, SIZE(molecule_kind_set)
873 71866 : molecule_kind => molecule_kind_set(i)
874 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
875 : molecule_list=molecule_list, &
876 : natom=natom, &
877 71866 : nbond=nbond, bond_list=bond_list)
878 71866 : molecule => molecule_set(molecule_list(1))
879 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
880 186817 : DO j = 1, nbond
881 114951 : atm_a = bond_list(j)%a
882 114951 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
883 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
884 114951 : name=name_atm_a)
885 114951 : atm_b = bond_list(j)%b
886 114951 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
887 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
888 114951 : name=name_atm_b)
889 114951 : found = .FALSE.
890 114951 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
891 114951 : CALL uppercase(name_atm_a)
892 114951 : CALL uppercase(name_atm_b)
893 :
894 : ! loop over params from GROMOS
895 114951 : 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 114951 : IF (ASSOCIATED(chm_info%bond_a)) THEN
913 1449364 : DO k = 1, SIZE(chm_info%bond_a)
914 : IF ((((chm_info%bond_a(k)) == (name_atm_a)) .AND. &
915 1449340 : ((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 41455 : bond_list(j)%bond_kind%id_type = do_ff_charmm
919 41455 : bond_list(j)%bond_kind%k(1) = chm_info%bond_k(k)
920 41455 : bond_list(j)%bond_kind%r0 = chm_info%bond_r0(k)
921 41455 : CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b)
922 41455 : found = .TRUE.
923 41455 : EXIT
924 : END IF
925 : END DO
926 : END IF
927 :
928 : ! loop over params from AMBER
929 114951 : 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 114951 : IF (ASSOCIATED(inp_info%bond_a)) THEN
947 10440 : DO k = 1, SIZE(inp_info%bond_a)
948 : IF ((((inp_info%bond_a(k)) == (name_atm_a)) .AND. &
949 10394 : ((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 8568 : bond_list(j)%bond_kind%id_type = inp_info%bond_kind(k)
953 59976 : bond_list(j)%bond_kind%k(:) = inp_info%bond_k(:, k)
954 8568 : bond_list(j)%bond_kind%r0 = inp_info%bond_r0(k)
955 8568 : bond_list(j)%bond_kind%cs = inp_info%bond_cs(k)
956 8568 : CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b)
957 8568 : found = .TRUE.
958 8568 : EXIT
959 : END IF
960 : END DO
961 : END IF
962 :
963 114951 : 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 186817 : 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 146369 : bond_list=bond_list)
977 :
978 : END DO
979 :
980 2637 : CALL timestop(handle2)
981 :
982 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1018 : LOGICAL :: found, only_qm
1019 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1020 2637 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
1021 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1022 : TYPE(molecule_type), POINTER :: molecule
1023 :
1024 2637 : CALL timeset(routineN, handle2)
1025 :
1026 2637 : IF (iw > 0) THEN
1027 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1028 242 : "FORCEFIELD| Checking for bend terms"
1029 : END IF
1030 :
1031 74503 : DO i = 1, SIZE(molecule_kind_set)
1032 71866 : molecule_kind => molecule_kind_set(i)
1033 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1034 : molecule_list=molecule_list, &
1035 : natom=natom, &
1036 71866 : nbend=nbend, bend_list=bend_list)
1037 71866 : molecule => molecule_set(molecule_list(1))
1038 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1039 212500 : DO j = 1, nbend
1040 140634 : atm_a = bend_list(j)%a
1041 140634 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1042 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1043 140634 : name=name_atm_a)
1044 140634 : atm_b = bend_list(j)%b
1045 140634 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1046 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1047 140634 : name=name_atm_b)
1048 140634 : atm_c = bend_list(j)%c
1049 140634 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1050 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1051 140634 : name=name_atm_c)
1052 140634 : found = .FALSE.
1053 140634 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c)
1054 140634 : CALL uppercase(name_atm_a)
1055 140634 : CALL uppercase(name_atm_b)
1056 140634 : CALL uppercase(name_atm_c)
1057 :
1058 : ! loop over params from GROMOS
1059 140634 : 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 140634 : IF (ASSOCIATED(chm_info%bend_a)) THEN
1076 6045175 : 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 6045101 : ((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 67527 : bend_list(j)%bend_kind%id_type = do_ff_charmm
1084 67527 : bend_list(j)%bend_kind%k = chm_info%bend_k(k)
1085 67527 : bend_list(j)%bend_kind%theta0 = chm_info%bend_theta0(k)
1086 : CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, &
1087 67527 : name_atm_c)
1088 67527 : found = .TRUE.
1089 67527 : EXIT
1090 : END IF
1091 : END DO
1092 : END IF
1093 :
1094 : ! loop over params from AMBER
1095 140634 : 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 140634 : IF (ASSOCIATED(inp_info%bend_a)) THEN
1116 28429 : 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 28413 : ((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 13411 : bend_list(j)%bend_kind%id_type = inp_info%bend_kind(k)
1124 13411 : bend_list(j)%bend_kind%k = inp_info%bend_k(k)
1125 13411 : bend_list(j)%bend_kind%theta0 = inp_info%bend_theta0(k)
1126 13411 : bend_list(j)%bend_kind%cb = inp_info%bend_cb(k)
1127 13411 : bend_list(j)%bend_kind%r012 = inp_info%bend_r012(k)
1128 13411 : bend_list(j)%bend_kind%r032 = inp_info%bend_r032(k)
1129 13411 : bend_list(j)%bend_kind%kbs12 = inp_info%bend_kbs12(k)
1130 13411 : bend_list(j)%bend_kind%kbs32 = inp_info%bend_kbs32(k)
1131 13411 : bend_list(j)%bend_kind%kss = inp_info%bend_kss(k)
1132 13411 : bend_list(j)%bend_kind%legendre%order = inp_info%bend_legendre(k)%order
1133 13411 : IF (bend_list(j)%bend_kind%legendre%order /= 0) THEN
1134 13411 : IF (ASSOCIATED(bend_list(j)%bend_kind%legendre%coeffs)) THEN
1135 11086 : DEALLOCATE (bend_list(j)%bend_kind%legendre%coeffs)
1136 : END IF
1137 40233 : ALLOCATE (bend_list(j)%bend_kind%legendre%coeffs(bend_list(j)%bend_kind%legendre%order))
1138 27062 : DO l = 1, bend_list(j)%bend_kind%legendre%order
1139 27062 : 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 13411 : name_atm_c)
1144 13411 : found = .TRUE.
1145 13411 : EXIT
1146 : END IF
1147 : END DO
1148 : END IF
1149 :
1150 140634 : 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 212500 : 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 146369 : bend_list=bend_list)
1164 : END DO
1165 2637 : CALL timestop(handle2)
1166 :
1167 2637 : 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 2637 : 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 2637 : 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 2637 : TYPE(ub_type), DIMENSION(:), POINTER :: ub_list
1202 :
1203 2637 : CALL timeset(routineN, handle2)
1204 :
1205 2637 : IF (iw > 0) THEN
1206 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1207 242 : "FORCEFIELD| Checking for Urey-Bradley (UB) terms"
1208 : END IF
1209 :
1210 74503 : DO i = 1, SIZE(molecule_kind_set)
1211 71866 : molecule_kind => molecule_kind_set(i)
1212 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1213 : molecule_list=molecule_list, &
1214 : natom=natom, &
1215 71866 : nub=nub, ub_list=ub_list)
1216 71866 : molecule => molecule_set(molecule_list(1))
1217 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1218 212342 : DO j = 1, nub
1219 140476 : atm_a = ub_list(j)%a
1220 140476 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1221 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1222 140476 : name=name_atm_a)
1223 140476 : atm_b = ub_list(j)%b
1224 140476 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1225 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1226 140476 : name=name_atm_b)
1227 140476 : atm_c = ub_list(j)%c
1228 140476 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1229 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1230 140476 : name=name_atm_c)
1231 140476 : found = .FALSE.
1232 140476 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c)
1233 140476 : CALL uppercase(name_atm_a)
1234 140476 : CALL uppercase(name_atm_b)
1235 140476 : 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 140476 : 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 140476 : IF (ASSOCIATED(inp_info%ub_a)) THEN
1270 50200 : 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 36781 : ((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 13419 : ((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 140476 : 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 119776 : array=Ainfo)
1299 119776 : ub_list(j)%id_type = do_ff_undef
1300 119776 : ub_list(j)%ub_kind%id_type = do_ff_undef
1301 479104 : ub_list(j)%ub_kind%k = 0.0_dp
1302 119776 : ub_list(j)%ub_kind%r0 = 0.0_dp
1303 : END IF
1304 :
1305 : ! QM/MM modifications
1306 212342 : 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 146369 : ub_list=ub_list)
1314 :
1315 : END DO
1316 :
1317 2637 : CALL timestop(handle2)
1318 :
1319 2637 : 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 2637 : 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 2637 : 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 2637 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
1362 :
1363 2637 : CALL timeset(routineN, handle2)
1364 :
1365 2637 : IF (iw > 0) THEN
1366 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1367 242 : "FORCEFIELD| Checking for torsion terms"
1368 : END IF
1369 :
1370 74503 : DO i = 1, SIZE(molecule_kind_set)
1371 71866 : 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 71866 : torsion_list=torsion_list)
1378 71866 : molecule => molecule_set(molecule_list(1))
1379 : CALL get_molecule(molecule=molecule, &
1380 : first_atom=first, &
1381 71866 : last_atom=last)
1382 229539 : DO j = 1, ntorsion
1383 229539 : IF (torsion_list(j)%torsion_kind%id_type == do_ff_undef) THEN
1384 116171 : atm_a = torsion_list(j)%a
1385 116171 : atomic_kind => particle_set(atm_a + first - 1)%atomic_kind
1386 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1387 116171 : name=name_atm_a)
1388 116171 : atm_b = torsion_list(j)%b
1389 116171 : atomic_kind => particle_set(atm_b + first - 1)%atomic_kind
1390 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1391 116171 : name=name_atm_b)
1392 116171 : atm_c = torsion_list(j)%c
1393 116171 : atomic_kind => particle_set(atm_c + first - 1)%atomic_kind
1394 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1395 116171 : name=name_atm_c)
1396 116171 : atm_d = torsion_list(j)%d
1397 116171 : atomic_kind => particle_set(atm_d + first - 1)%atomic_kind
1398 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
1399 116171 : name=name_atm_d)
1400 116171 : found = .FALSE.
1401 116171 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d)
1402 116171 : CALL uppercase(name_atm_a)
1403 116171 : CALL uppercase(name_atm_b)
1404 116171 : CALL uppercase(name_atm_c)
1405 116171 : CALL uppercase(name_atm_d)
1406 :
1407 : ! Loop over params from GROMOS
1408 116171 : 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 116171 : 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 116171 : 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 116171 : 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 116171 : 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 36082 : array=Ainfo)
1576 36082 : torsion_list(j)%torsion_kind%id_type = do_ff_undef
1577 36082 : torsion_list(j)%id_type = do_ff_undef
1578 : END IF
1579 :
1580 : ! QM/MM modifications
1581 116171 : 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 146369 : torsion_list=torsion_list)
1601 :
1602 : END DO ! molecule kind
1603 :
1604 2637 : CALL timestop(handle2)
1605 :
1606 2637 : 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 2637 : 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 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list
1640 : LOGICAL :: found, only_qm
1641 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1642 2637 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list
1643 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1644 : TYPE(molecule_type), POINTER :: molecule
1645 :
1646 2637 : CALL timeset(routineN, handle2)
1647 :
1648 2637 : IF (iw > 0) THEN
1649 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1650 242 : "FORCEFIELD| Checking for improper terms"
1651 : END IF
1652 :
1653 74503 : DO i = 1, SIZE(molecule_kind_set)
1654 :
1655 71866 : 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 71866 : impr_list=impr_list)
1661 :
1662 71866 : molecule => molecule_set(molecule_list(1))
1663 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1664 :
1665 77178 : 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 77178 : 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 146369 : CALL set_molecule_kind(molecule_kind=molecule_kind, impr_list=impr_list)
1814 :
1815 : END DO
1816 :
1817 2637 : CALL timestop(handle2)
1818 :
1819 2637 : 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 2637 : 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 2637 : 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 2637 : TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list
1856 :
1857 2637 : CALL timeset(routineN, handle2)
1858 :
1859 2637 : IF (iw > 0) THEN
1860 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1861 242 : "FORCEFIELD| Checking for out-of-plane bend terms"
1862 : END IF
1863 :
1864 74503 : DO i = 1, SIZE(molecule_kind_set)
1865 71866 : molecule_kind => molecule_kind_set(i)
1866 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1867 : molecule_list=molecule_list, &
1868 : natom=natom, &
1869 71866 : nopbend=nopbend, opbend_list=opbend_list)
1870 71866 : molecule => molecule_set(molecule_list(1))
1871 :
1872 71866 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
1873 77178 : 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 77178 : 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 146369 : CALL set_molecule_kind(molecule_kind=molecule_kind, opbend_list=opbend_list)
1953 :
1954 : END DO
1955 :
1956 2637 : CALL timestop(handle2)
1957 :
1958 2637 : 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 : 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)) THEN
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 : END IF
2035 36 : charge = -HUGE(0.0_dp)
2036 :
2037 : ! Check if the potential really requires the charge definition..
2038 36 : IF (ASSOCIATED(inp_info%nonbonded)) THEN
2039 18 : IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN
2040 : ! Let's find the nonbonded potential where this atom is involved
2041 18 : only_manybody = .TRUE.
2042 18 : found_p = .FALSE.
2043 30 : DO j = 1, SIZE(inp_info%nonbonded%pot)
2044 30 : IF (atmname == inp_info%nonbonded%pot(j)%pot%at1 .OR. &
2045 0 : atmname == inp_info%nonbonded%pot(j)%pot%at2) THEN
2046 18 : SELECT CASE (inp_info%nonbonded%pot(j)%pot%type(1))
2047 : CASE (ea_type, tersoff_type, siepmann_type)
2048 : ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential
2049 : ! Do nothing..
2050 : CASE DEFAULT
2051 : only_manybody = .FALSE.
2052 18 : EXIT
2053 : END SELECT
2054 : found_p = .TRUE.
2055 : END IF
2056 : END DO
2057 18 : IF (only_manybody .AND. found_p) THEN
2058 0 : charges(iatom) = 0.0_dp
2059 : END IF
2060 : END IF
2061 : END IF
2062 :
2063 : ! QM/MM modifications
2064 80 : IF (only_qm .AND. my_qmmm) THEN
2065 6 : IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN
2066 6 : scale_factor = 0.0_dp
2067 6 : IF (is_link_atom) THEN
2068 : ! Find the scaling factor...
2069 0 : DO ilink = 1, SIZE(qmmm_env%mm_link_atoms)
2070 0 : IF (iatom == qmmm_env%mm_link_atoms(ilink)) EXIT
2071 : END DO
2072 0 : CPASSERT(ilink <= SIZE(qmmm_env%mm_link_atoms))
2073 0 : scale_factor = qmmm_env%fist_scale_charge_link(ilink)
2074 : END IF
2075 6 : charges(iatom) = charges(iatom)*scale_factor
2076 : END IF
2077 : END IF
2078 : END DO
2079 :
2080 : ! Sum up total charges for IO
2081 44 : charge_tot = SUM(charges)
2082 :
2083 : ! Print total charge of the system
2084 8 : IF (iw4 > 0) THEN
2085 : WRITE (UNIT=iw4, FMT="(/,T2,A,T61,F20.10)") &
2086 4 : "FORCEFIELD| Total charge of the classical system: ", charge_tot
2087 : END IF
2088 :
2089 8 : CALL timestop(handle)
2090 :
2091 16 : END SUBROUTINE force_field_pack_charges
2092 :
2093 : ! **************************************************************************************************
2094 : !> \brief Set up atomic_kind_set()%fist_potential%[qeff]
2095 : !> and shell potential parameters
2096 : !> \param atomic_kind_set ...
2097 : !> \param qmmm_env ...
2098 : !> \param fatal ...
2099 : !> \param iw ...
2100 : !> \param iw4 ...
2101 : !> \param Ainfo ...
2102 : !> \param my_qmmm ...
2103 : !> \param inp_info ...
2104 : ! **************************************************************************************************
2105 2629 : SUBROUTINE force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4, &
2106 : Ainfo, my_qmmm, inp_info)
2107 :
2108 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2109 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
2110 : LOGICAL, INTENT(INOUT) :: fatal
2111 : INTEGER, INTENT(IN) :: iw, iw4
2112 : CHARACTER(LEN=default_string_length), &
2113 : DIMENSION(:), POINTER :: Ainfo
2114 : LOGICAL, INTENT(IN) :: my_qmmm
2115 : TYPE(input_info_type), POINTER :: inp_info
2116 :
2117 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charge'
2118 :
2119 : CHARACTER(LEN=default_string_length) :: atmname
2120 : INTEGER :: handle, i, ilink, j
2121 2629 : INTEGER, DIMENSION(:), POINTER :: my_atom_list
2122 : LOGICAL :: found, found_p, is_link_atom, is_shell, &
2123 : only_manybody, only_qm
2124 : REAL(KIND=dp) :: charge, charge_tot, cs_charge, &
2125 : scale_factor
2126 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2127 : TYPE(fist_potential_type), POINTER :: fist_potential
2128 :
2129 2629 : CALL timeset(routineN, handle)
2130 :
2131 2629 : charge_tot = 0.0_dp
2132 :
2133 13872 : DO i = 1, SIZE(atomic_kind_set)
2134 11243 : atomic_kind => atomic_kind_set(i)
2135 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2136 : fist_potential=fist_potential, &
2137 : atom_list=my_atom_list, &
2138 11243 : name=atmname)
2139 11243 : CALL get_potential(potential=fist_potential, qeff=charge)
2140 :
2141 11243 : is_shell = .FALSE.
2142 11243 : found = .FALSE.
2143 11243 : only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom)
2144 11243 : CALL uppercase(atmname)
2145 11243 : IF (charge /= -HUGE(0.0_dp)) found = .TRUE.
2146 :
2147 : ! Always have the input param last to overwrite all the other ones
2148 11243 : IF (ASSOCIATED(inp_info%charge_atm)) THEN
2149 5559 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") ""
2150 27280 : DO j = 1, SIZE(inp_info%charge_atm)
2151 : IF (debug_this_module) THEN
2152 : IF (iw > 0) THEN
2153 : WRITE (UNIT=iw, FMT="(T2,A)") &
2154 : "Checking charges for the atomic kinds "// &
2155 : TRIM(inp_info%charge_atm(j))//" and "//TRIM(atmname)
2156 : END IF
2157 : END IF
2158 27280 : IF ((inp_info%charge_atm(j)) == atmname) THEN
2159 5471 : charge = inp_info%charge(j)
2160 5471 : CALL issue_duplications(found, "Charge", atmname)
2161 5471 : found = .TRUE.
2162 : END IF
2163 : END DO
2164 : END IF
2165 : ! Check if the ATOM type has a core-shell associated.. In this case
2166 : ! print a warning: the CHARGE will not be used if defined.. or we can
2167 : ! even skip its definition..
2168 11243 : IF (ASSOCIATED(inp_info%shell_list)) THEN
2169 1410 : DO j = 1, SIZE(inp_info%shell_list)
2170 1410 : IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN
2171 450 : is_shell = .TRUE.
2172 : cs_charge = inp_info%shell_list(j)%shell%charge_core + &
2173 450 : inp_info%shell_list(j)%shell%charge_shell
2174 450 : charge = 0.0_dp
2175 450 : IF (found) THEN
2176 : IF (found) THEN
2177 : CALL cp_warn(__LOCATION__, &
2178 : "CORE-SHELL model defined for KIND ("//TRIM(atmname)//")"// &
2179 204 : " ignoring charge definition! ")
2180 : END IF
2181 : ELSE
2182 246 : found = .TRUE.
2183 : END IF
2184 : END IF
2185 : END DO
2186 : END IF
2187 : ! Check if the potential really requires the charge definition..
2188 11243 : IF (ASSOCIATED(inp_info%nonbonded)) THEN
2189 4323 : IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN
2190 : ! Let's find the nonbonded potential where this atom is involved
2191 4323 : only_manybody = .TRUE.
2192 4323 : found_p = .FALSE.
2193 7841 : DO j = 1, SIZE(inp_info%nonbonded%pot)
2194 7644 : IF (atmname == inp_info%nonbonded%pot(j)%pot%at1 .OR. &
2195 197 : atmname == inp_info%nonbonded%pot(j)%pot%at2) THEN
2196 4312 : SELECT CASE (inp_info%nonbonded%pot(j)%pot%type(1))
2197 : CASE (ea_type, tersoff_type, siepmann_type, nequip_type, &
2198 : allegro_type, deepmd_type, ace_type, mace_type)
2199 : ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential
2200 : ! Do nothing..
2201 : CASE DEFAULT
2202 : only_manybody = .FALSE.
2203 4312 : EXIT
2204 : END SELECT
2205 : found_p = .TRUE.
2206 : END IF
2207 : END DO
2208 4323 : IF (only_manybody .AND. found_p) THEN
2209 150 : charge = 0.0_dp
2210 150 : found = .TRUE.
2211 : END IF
2212 : END IF
2213 : END IF
2214 11243 : IF (.NOT. found) THEN
2215 : ! Set the charge to zero anyway in case the user decides to ignore
2216 : ! missing critical parameters.
2217 12 : charge = 0.0_dp
2218 : CALL store_FF_missing_par(atm1=TRIM(atmname), &
2219 : fatal=fatal, &
2220 : type_name="Charge", &
2221 12 : array=Ainfo)
2222 : END IF
2223 : !
2224 : ! QM/MM modifications
2225 : !
2226 11243 : IF (only_qm .AND. my_qmmm) THEN
2227 1286 : IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN
2228 1076 : scale_factor = 0.0_dp
2229 1076 : IF (is_link_atom) THEN
2230 : !
2231 : ! Find the scaling factor...
2232 : !
2233 386 : DO ilink = 1, SIZE(qmmm_env%mm_link_atoms)
2234 658 : IF (ANY(my_atom_list == qmmm_env%mm_link_atoms(ilink))) EXIT
2235 : END DO
2236 114 : CPASSERT(ilink <= SIZE(qmmm_env%mm_link_atoms))
2237 114 : scale_factor = qmmm_env%fist_scale_charge_link(ilink)
2238 : END IF
2239 1076 : charge = charge*scale_factor
2240 : END IF
2241 : END IF
2242 :
2243 11243 : CALL set_potential(potential=fist_potential, qeff=charge)
2244 : ! Sum up total charges for IO
2245 13872 : IF (found) THEN
2246 11231 : IF (is_shell) THEN
2247 450 : charge_tot = charge_tot + atomic_kind%natom*cs_charge
2248 : ELSE
2249 10781 : charge_tot = charge_tot + atomic_kind%natom*charge
2250 : END IF
2251 : END IF
2252 : END DO
2253 :
2254 : ! Print total charge of the system
2255 2629 : IF (iw4 > 0) THEN
2256 : WRITE (UNIT=iw4, FMT="(/,T2,A,T61,F20.10)") &
2257 1300 : "FORCEFIELD| Total charge of the classical system: ", charge_tot
2258 : END IF
2259 :
2260 2629 : CALL timestop(handle)
2261 :
2262 2629 : END SUBROUTINE force_field_pack_charge
2263 :
2264 : ! **************************************************************************************************
2265 : !> \brief Set up the radius of the electrostatic multipole in Fist
2266 : !> \param atomic_kind_set ...
2267 : !> \param iw ...
2268 : !> \param subsys_section ...
2269 : !> \author Toon.Verstraelen@gmail.com
2270 : ! **************************************************************************************************
2271 5274 : SUBROUTINE force_field_pack_radius(atomic_kind_set, iw, subsys_section)
2272 :
2273 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2274 : INTEGER, INTENT(IN) :: iw
2275 : TYPE(section_vals_type), POINTER :: subsys_section
2276 :
2277 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_radius'
2278 :
2279 : CHARACTER(LEN=default_string_length) :: inp_kind_name, kind_name
2280 : INTEGER :: handle, i, i_rep, n_rep
2281 : LOGICAL :: found
2282 : REAL(KIND=dp) :: mm_radius
2283 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2284 : TYPE(fist_potential_type), POINTER :: fist_potential
2285 : TYPE(section_vals_type), POINTER :: kind_section
2286 :
2287 2637 : CALL timeset(routineN, handle)
2288 :
2289 2637 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
2290 2637 : CALL section_vals_get(kind_section, n_repetition=n_rep)
2291 :
2292 13900 : DO i = 1, SIZE(atomic_kind_set)
2293 11263 : atomic_kind => atomic_kind_set(i)
2294 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2295 11263 : fist_potential=fist_potential, name=kind_name)
2296 11263 : CALL uppercase(kind_name)
2297 11263 : found = .FALSE.
2298 :
2299 : ! Try to find a matching KIND section in the SUBSYS section and read the
2300 : ! MM_RADIUS field if it is present. In case the kind section is never
2301 : ! encountered, the mm_radius remains zero.
2302 11263 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") ""
2303 11263 : mm_radius = 0.0_dp
2304 39554 : DO i_rep = 1, n_rep
2305 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
2306 28291 : c_val=inp_kind_name, i_rep_section=i_rep)
2307 28291 : CALL uppercase(inp_kind_name)
2308 28291 : IF (iw > 0) THEN
2309 : WRITE (UNIT=iw, FMT="(T2,A)") &
2310 : "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// &
2311 905 : " and "//TRIM(inp_kind_name)//" for MM_RADIUS"
2312 : END IF
2313 39554 : IF (TRIM(kind_name) == TRIM(inp_kind_name)) THEN
2314 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
2315 1839 : keyword_name="MM_RADIUS", r_val=mm_radius)
2316 1839 : CALL issue_duplications(found, "MM_RADIUS", kind_name)
2317 1839 : found = .TRUE.
2318 : END IF
2319 : END DO
2320 13900 : CALL set_potential(potential=fist_potential, mm_radius=mm_radius)
2321 : END DO
2322 :
2323 2637 : CALL timestop(handle)
2324 :
2325 2637 : END SUBROUTINE force_field_pack_radius
2326 :
2327 : ! **************************************************************************************************
2328 : !> \brief Set up the polarizable FF parameters
2329 : !> \param atomic_kind_set ...
2330 : !> \param iw ...
2331 : !> \param inp_info ...
2332 : !> \author Toon.Verstraelen@gmail.com
2333 : ! **************************************************************************************************
2334 2637 : SUBROUTINE force_field_pack_pol(atomic_kind_set, iw, inp_info)
2335 :
2336 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2337 : INTEGER, INTENT(IN) :: iw
2338 : TYPE(input_info_type), POINTER :: inp_info
2339 :
2340 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_pol'
2341 :
2342 : CHARACTER(LEN=default_string_length) :: kind_name
2343 : INTEGER :: handle, i, j
2344 : LOGICAL :: found
2345 : REAL(KIND=dp) :: apol, cpol
2346 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2347 : TYPE(fist_potential_type), POINTER :: fist_potential
2348 :
2349 2637 : CALL timeset(routineN, handle)
2350 :
2351 2637 : IF (iw > 0) THEN
2352 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2353 242 : "FORCEFIELD| Checking for polarisable forcefield terms"
2354 : END IF
2355 :
2356 13900 : DO i = 1, SIZE(atomic_kind_set)
2357 11263 : atomic_kind => atomic_kind_set(i)
2358 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2359 : fist_potential=fist_potential, &
2360 11263 : name=kind_name)
2361 11263 : CALL get_potential(potential=fist_potential, apol=apol, cpol=cpol)
2362 11263 : CALL uppercase(kind_name)
2363 11263 : found = .FALSE.
2364 :
2365 11263 : IF (iw > 0) WRITE (UNIT=iw, FMT="(A)") ""
2366 : ! Always have the input param last to overwrite all the other ones
2367 11263 : IF (ASSOCIATED(inp_info%apol_atm)) THEN
2368 292 : DO j = 1, SIZE(inp_info%apol_atm)
2369 200 : IF (iw > 0) THEN
2370 : WRITE (UNIT=iw, FMT="(T2,A)") &
2371 : "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// &
2372 0 : " and "//TRIM(inp_info%apol_atm(j))//" for APOL"
2373 : END IF
2374 292 : IF ((inp_info%apol_atm(j)) == kind_name) THEN
2375 64 : apol = inp_info%apol(j)
2376 64 : CALL issue_duplications(found, "APOL", kind_name)
2377 64 : found = .TRUE.
2378 : END IF
2379 : END DO
2380 : END IF
2381 :
2382 11263 : IF (ASSOCIATED(inp_info%cpol_atm)) THEN
2383 0 : DO j = 1, SIZE(inp_info%cpol_atm)
2384 0 : IF (iw > 0) THEN
2385 : WRITE (UNIT=iw, FMT="(T2,A)") &
2386 : "FORCEFIELD| Matching atomic kinds "//TRIM(kind_name)// &
2387 0 : " and "//TRIM(inp_info%cpol_atm(j))//" for CPOL"
2388 : END IF
2389 0 : IF ((inp_info%cpol_atm(j)) == kind_name) THEN
2390 0 : cpol = inp_info%cpol(j)
2391 0 : CALL issue_duplications(found, "CPOL", kind_name)
2392 0 : found = .TRUE.
2393 : END IF
2394 : END DO
2395 : END IF
2396 :
2397 13900 : CALL set_potential(potential=fist_potential, apol=apol, cpol=cpol)
2398 :
2399 : END DO
2400 :
2401 2637 : CALL timestop(handle)
2402 :
2403 2637 : END SUBROUTINE force_field_pack_pol
2404 :
2405 : ! **************************************************************************************************
2406 : !> \brief Set up damping parameters
2407 : !> \param atomic_kind_set ...
2408 : !> \param iw ...
2409 : !> \param inp_info ...
2410 : ! **************************************************************************************************
2411 2637 : SUBROUTINE force_field_pack_damp(atomic_kind_set, iw, inp_info)
2412 :
2413 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2414 : INTEGER :: iw
2415 : TYPE(input_info_type), POINTER :: inp_info
2416 :
2417 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_damp'
2418 :
2419 : CHARACTER(len=default_string_length) :: atm_name1, atm_name2, my_atm_name1, &
2420 : my_atm_name2
2421 : INTEGER :: handle2, i, j, k, nkinds
2422 : LOGICAL :: found
2423 : TYPE(atomic_kind_type), POINTER :: atomic_kind, atomic_kind2
2424 : TYPE(damping_p_type), POINTER :: damping
2425 :
2426 2637 : CALL timeset(routineN, handle2)
2427 :
2428 2637 : IF (iw > 0) THEN
2429 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2430 242 : "FORCEFIELD| Checking for damping terms"
2431 : END IF
2432 :
2433 2637 : NULLIFY (damping)
2434 2637 : nkinds = SIZE(atomic_kind_set)
2435 :
2436 13900 : DO j = 1, SIZE(atomic_kind_set)
2437 :
2438 11263 : atomic_kind => atomic_kind_set(j)
2439 :
2440 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2441 11263 : name=atm_name1)
2442 11263 : CALL uppercase(atm_name1)
2443 :
2444 11263 : IF (ASSOCIATED(inp_info%damping_list)) THEN
2445 50 : DO i = 1, SIZE(inp_info%damping_list)
2446 28 : my_atm_name1 = inp_info%damping_list(i)%atm_name1
2447 28 : my_atm_name2 = inp_info%damping_list(i)%atm_name2
2448 : IF (debug_this_module) THEN
2449 : IF (iw > 0) THEN
2450 : WRITE (UNIT=iw, FMT="(T2,A)") &
2451 : "FORCEFIELD| Check damping for the atomic kinds "// &
2452 : TRIM(my_atm_name1)//" and "//TRIM(atm_name1)
2453 : END IF
2454 : END IF
2455 50 : IF (my_atm_name1 == atm_name1) THEN
2456 12 : IF (.NOT. ASSOCIATED(damping)) THEN
2457 10 : CALL damping_p_create(damping, nkinds)
2458 : END IF
2459 12 : found = .FALSE.
2460 40 : DO k = 1, SIZE(atomic_kind_set)
2461 28 : atomic_kind2 => atomic_kind_set(k)
2462 : CALL get_atomic_kind(atomic_kind=atomic_kind2, &
2463 28 : name=atm_name2)
2464 28 : CALL uppercase(atm_name2)
2465 40 : IF (my_atm_name2 == atm_name2) THEN
2466 12 : IF (damping%damp(k)%bij /= HUGE(0.0_dp)) found = .TRUE.
2467 12 : CALL issue_duplications(found, "Damping", atm_name1)
2468 12 : found = .TRUE.
2469 24 : SELECT CASE (TRIM(inp_info%damping_list(i)%dtype))
2470 : CASE ('TANG-TOENNIES')
2471 12 : damping%damp(k)%itype = tang_toennies
2472 : CASE DEFAULT
2473 24 : CPABORT("Unknown damping type.")
2474 : END SELECT
2475 12 : damping%damp(k)%order = inp_info%damping_list(i)%order
2476 12 : damping%damp(k)%bij = inp_info%damping_list(i)%bij
2477 12 : damping%damp(k)%cij = inp_info%damping_list(i)%cij
2478 : END IF
2479 : END DO
2480 12 : IF (.NOT. found) THEN
2481 : CALL cp_warn(__LOCATION__, &
2482 : "Atom "//TRIM(my_atm_name2)// &
2483 : " in damping parameters for atom "//TRIM(my_atm_name1)// &
2484 0 : " not found.")
2485 : END IF
2486 : END IF
2487 : END DO
2488 : END IF
2489 :
2490 11263 : CALL set_atomic_kind(atomic_kind=atomic_kind, damping=damping)
2491 :
2492 13900 : NULLIFY (damping)
2493 :
2494 : END DO
2495 :
2496 2637 : CALL timestop(handle2)
2497 :
2498 2637 : END SUBROUTINE force_field_pack_damp
2499 :
2500 : ! **************************************************************************************************
2501 : !> \brief Set up shell potential parameters
2502 : !> \param particle_set ...
2503 : !> \param atomic_kind_set ...
2504 : !> \param molecule_kind_set ...
2505 : !> \param molecule_set ...
2506 : !> \param root_section ...
2507 : !> \param subsys_section ...
2508 : !> \param shell_particle_set ...
2509 : !> \param core_particle_set ...
2510 : !> \param cell ...
2511 : !> \param iw ...
2512 : !> \param inp_info ...
2513 : ! **************************************************************************************************
2514 13185 : SUBROUTINE force_field_pack_shell(particle_set, atomic_kind_set, &
2515 : molecule_kind_set, molecule_set, root_section, subsys_section, &
2516 : shell_particle_set, core_particle_set, cell, iw, inp_info)
2517 :
2518 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2519 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2520 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2521 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
2522 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
2523 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_particle_set, core_particle_set
2524 : TYPE(cell_type), POINTER :: cell
2525 : INTEGER :: iw
2526 : TYPE(input_info_type), POINTER :: inp_info
2527 :
2528 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_shell'
2529 :
2530 : CHARACTER(LEN=default_string_length) :: atmname
2531 : INTEGER :: counter, first, first_shell, handle2, i, &
2532 : j, last, last_shell, n, natom, nmol, &
2533 : nshell_tot
2534 2637 : INTEGER, DIMENSION(:), POINTER :: molecule_list, shell_list_tmp
2535 : LOGICAL :: core_coord_read, found_shell, is_a_shell, is_link_atom, null_massfrac, only_qm, &
2536 : save_mem, shell_adiabatic, shell_coord_read
2537 : REAL(KIND=dp) :: atmmass
2538 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2539 : TYPE(molecule_kind_type), POINTER :: molecule_kind
2540 : TYPE(molecule_type), POINTER :: molecule
2541 : TYPE(section_vals_type), POINTER :: global_section
2542 : TYPE(shell_kind_type), POINTER :: shell
2543 2637 : TYPE(shell_type), DIMENSION(:), POINTER :: shell_list
2544 :
2545 2637 : CALL timeset(routineN, handle2)
2546 :
2547 2637 : nshell_tot = 0
2548 2637 : n = 0
2549 2637 : first_shell = 1
2550 2637 : null_massfrac = .FALSE.
2551 2637 : core_coord_read = .FALSE.
2552 2637 : shell_coord_read = .FALSE.
2553 :
2554 2637 : NULLIFY (global_section)
2555 2637 : global_section => section_vals_get_subs_vals(root_section, "GLOBAL")
2556 2637 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
2557 :
2558 2637 : IF (iw > 0) THEN
2559 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2560 242 : "FORCEFIELD| Checking for core-shell terms"
2561 : END IF
2562 :
2563 13900 : DO i = 1, SIZE(atomic_kind_set)
2564 11263 : atomic_kind => atomic_kind_set(i)
2565 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2566 11263 : name=atmname)
2567 :
2568 11263 : found_shell = .FALSE.
2569 11263 : only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom)
2570 11263 : CALL uppercase(atmname)
2571 :
2572 : ! The shell potential can be defined only from input
2573 13900 : IF (ASSOCIATED(inp_info%shell_list)) THEN
2574 1410 : DO j = 1, SIZE(inp_info%shell_list)
2575 : IF (debug_this_module) THEN
2576 : IF (iw > 0) THEN
2577 : WRITE (UNIT=iw, FMT="(T2,A)") &
2578 : "Checking shells for the atomic kinds "// &
2579 : TRIM(inp_info%shell_list(j)%atm_name)//" and "//TRIM(atmname)
2580 : END IF
2581 : END IF
2582 1410 : IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN
2583 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
2584 450 : shell=shell, mass=atmmass, natom=natom)
2585 450 : IF (.NOT. ASSOCIATED(shell)) ALLOCATE (shell)
2586 450 : nshell_tot = nshell_tot + natom
2587 450 : shell%charge_core = inp_info%shell_list(j)%shell%charge_core
2588 450 : shell%charge_shell = inp_info%shell_list(j)%shell%charge_shell
2589 450 : shell%massfrac = inp_info%shell_list(j)%shell%massfrac
2590 450 : IF (shell%massfrac < EPSILON(1.0_dp)) null_massfrac = .TRUE.
2591 450 : shell%k2_spring = inp_info%shell_list(j)%shell%k2_spring
2592 450 : shell%k4_spring = inp_info%shell_list(j)%shell%k4_spring
2593 450 : shell%max_dist = inp_info%shell_list(j)%shell%max_dist
2594 450 : shell%shell_cutoff = inp_info%shell_list(j)%shell%shell_cutoff
2595 450 : shell%mass_shell = shell%massfrac*atmmass
2596 450 : shell%mass_core = atmmass - shell%mass_shell
2597 450 : CALL issue_duplications(found_shell, "Shell", atmname)
2598 450 : found_shell = .TRUE.
2599 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
2600 450 : shell=shell, shell_active=.TRUE.)
2601 : END IF
2602 : END DO ! shell kind
2603 : END IF ! associated shell_list
2604 : END DO ! atomic kind
2605 :
2606 2637 : IF (iw > 0) THEN
2607 : WRITE (UNIT=iw, FMT="(/,T2,A,T61,I20)") &
2608 242 : "FORCEFIELD| Total number of particles with a shell:", nshell_tot
2609 : END IF
2610 : ! If shell-model is present: Create particle_set of shells (coord. vel. force)
2611 2637 : NULLIFY (shell_particle_set)
2612 2637 : NULLIFY (core_particle_set)
2613 2637 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, shell_adiabatic=shell_adiabatic)
2614 2637 : IF (nshell_tot > 0) THEN
2615 252 : IF (shell_adiabatic .AND. null_massfrac) THEN
2616 0 : CPABORT("Shell-model adiabatic: at least one shell_kind has mass zero")
2617 : END IF
2618 252 : CALL allocate_particle_set(shell_particle_set, nshell_tot)
2619 252 : CALL allocate_particle_set(core_particle_set, nshell_tot)
2620 252 : counter = 0
2621 : ! Initialise the shell (and core) coordinates with the particle (atomic) coordinates,
2622 : ! count the shell and set pointers
2623 28702 : DO i = 1, SIZE(particle_set)
2624 28450 : NULLIFY (atomic_kind)
2625 28450 : NULLIFY (shell)
2626 28450 : atomic_kind => particle_set(i)%atomic_kind
2627 28450 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_a_shell)
2628 28702 : IF (is_a_shell) THEN
2629 27962 : counter = counter + 1
2630 27962 : particle_set(i)%shell_index = counter
2631 27962 : shell_particle_set(counter)%shell_index = counter
2632 27962 : shell_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind
2633 195734 : shell_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3)
2634 27962 : shell_particle_set(counter)%atom_index = i
2635 27962 : core_particle_set(counter)%shell_index = counter
2636 27962 : core_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind
2637 195734 : core_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3)
2638 27962 : core_particle_set(counter)%atom_index = i
2639 : ELSE
2640 488 : particle_set(i)%shell_index = 0
2641 : END IF
2642 : END DO
2643 252 : CPASSERT(counter == nshell_tot)
2644 : END IF
2645 :
2646 : ! Read the shell (and core) coordinates from the restart file, if available
2647 : CALL read_binary_cs_coordinates("SHELL", shell_particle_set, root_section, &
2648 2637 : subsys_section, shell_coord_read, cell)
2649 : CALL read_binary_cs_coordinates("CORE", core_particle_set, root_section, &
2650 2637 : subsys_section, core_coord_read, cell)
2651 :
2652 2637 : IF (nshell_tot > 0) THEN
2653 : ! Read the shell (and core) coordinates from the input, if no coordinates were found
2654 : ! in the restart file
2655 252 : IF (shell_adiabatic) THEN
2656 252 : IF (.NOT. (core_coord_read .AND. shell_coord_read)) THEN
2657 : CALL read_shell_coord_input(particle_set, shell_particle_set, cell, &
2658 : subsys_section, core_particle_set, &
2659 236 : save_mem=save_mem)
2660 : END IF
2661 : ELSE
2662 0 : IF (.NOT. shell_coord_read) THEN
2663 : CALL read_shell_coord_input(particle_set, shell_particle_set, cell, &
2664 0 : subsys_section, save_mem=save_mem)
2665 : END IF
2666 : END IF
2667 : ! Determine the number of shells per molecule kind
2668 252 : n = 0
2669 11316 : DO i = 1, SIZE(molecule_kind_set)
2670 11064 : molecule_kind => molecule_kind_set(i)
2671 : CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, &
2672 11064 : natom=natom, nmolecule=nmol)
2673 11064 : molecule => molecule_set(molecule_list(1))
2674 11064 : CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last)
2675 33192 : ALLOCATE (shell_list_tmp(natom))
2676 11064 : counter = 0
2677 23106 : DO j = first, last
2678 12042 : atomic_kind => particle_set(j)%atomic_kind
2679 12042 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_a_shell)
2680 23106 : IF (is_a_shell) THEN
2681 11696 : counter = counter + 1
2682 11696 : shell_list_tmp(counter) = j - first + 1
2683 11696 : first_shell = MIN(first_shell, MAX(1, particle_set(j)%shell_index))
2684 : END IF
2685 : END DO ! j atom in molecule_kind i, molecule 1 of the molecule_list
2686 11064 : IF (counter /= 0) THEN
2687 : ! Setup of fist_shell and last_shell for all molecules..
2688 29072 : DO j = 1, SIZE(molecule_list)
2689 18306 : last_shell = first_shell + counter - 1
2690 18306 : molecule => molecule_set(molecule_list(j))
2691 18306 : molecule%first_shell = first_shell
2692 18306 : molecule%last_shell = last_shell
2693 29072 : first_shell = last_shell + 1
2694 : END DO
2695 : ! Setup of shell_list
2696 10766 : CALL get_molecule_kind(molecule_kind=molecule_kind, shell_list=shell_list)
2697 10766 : IF (ASSOCIATED(shell_list)) THEN
2698 0 : DEALLOCATE (shell_list)
2699 : END IF
2700 43994 : ALLOCATE (shell_list(counter))
2701 22462 : DO j = 1, counter
2702 11696 : shell_list(j)%a = shell_list_tmp(j)
2703 11696 : atomic_kind => particle_set(shell_list_tmp(j) + first - 1)%atomic_kind
2704 11696 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname, shell=shell)
2705 11696 : CALL uppercase(atmname)
2706 11696 : shell_list(j)%name = atmname
2707 22462 : shell_list(j)%shell_kind => shell
2708 : END DO
2709 10766 : CALL set_molecule_kind(molecule_kind=molecule_kind, nshell=counter, shell_list=shell_list)
2710 : END IF
2711 11064 : DEALLOCATE (shell_list_tmp)
2712 22380 : n = n + nmol*counter
2713 : END DO ! i molecule kind
2714 : END IF
2715 :
2716 2637 : CPASSERT(first_shell - 1 == nshell_tot)
2717 2637 : CPASSERT(n == nshell_tot)
2718 :
2719 2637 : CALL timestop(handle2)
2720 :
2721 2637 : END SUBROUTINE force_field_pack_shell
2722 :
2723 : ! **************************************************************************************************
2724 : !> \brief Assign input and potential info to potparm_nonbond14
2725 : !> \param atomic_kind_set ...
2726 : !> \param ff_type ...
2727 : !> \param qmmm_env ...
2728 : !> \param iw ...
2729 : !> \param Ainfo ...
2730 : !> \param chm_info ...
2731 : !> \param inp_info ...
2732 : !> \param gro_info ...
2733 : !> \param amb_info ...
2734 : !> \param potparm_nonbond14 ...
2735 : !> \param ewald_env ...
2736 : ! **************************************************************************************************
2737 2621 : SUBROUTINE force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, &
2738 : Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env)
2739 :
2740 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2741 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
2742 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
2743 : INTEGER :: iw
2744 : CHARACTER(LEN=default_string_length), &
2745 : DIMENSION(:), POINTER :: Ainfo
2746 : TYPE(charmm_info_type), POINTER :: chm_info
2747 : TYPE(input_info_type), POINTER :: inp_info
2748 : TYPE(gromos_info_type), POINTER :: gro_info
2749 : TYPE(amber_info_type), POINTER :: amb_info
2750 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond14
2751 : TYPE(ewald_environment_type), POINTER :: ewald_env
2752 :
2753 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond14'
2754 :
2755 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
2756 : name_atm_b, name_atm_b_local
2757 : INTEGER :: handle2, i, ii, j, jj, k, match_names
2758 : LOGICAL :: found, found_a, found_b, only_qm, &
2759 : use_qmmm_ff
2760 : REAL(KIND=dp) :: epsilon0, epsilon_a, epsilon_b, &
2761 : ewald_rcut, rmin, rmin2_a, rmin2_b
2762 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2763 : TYPE(pair_potential_single_type), POINTER :: pot
2764 :
2765 2621 : CALL timeset(routineN, handle2)
2766 :
2767 2621 : use_qmmm_ff = qmmm_env%use_qmmm_ff
2768 2621 : NULLIFY (pot)
2769 :
2770 2621 : IF (iw > 0) THEN
2771 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
2772 241 : "FORCEFIELD| Checking for nonbonded14 terms"
2773 : END IF
2774 :
2775 2621 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
2776 2621 : CALL pair_potential_pp_create(potparm_nonbond14, SIZE(atomic_kind_set))
2777 :
2778 13810 : DO i = 1, SIZE(atomic_kind_set)
2779 11189 : atomic_kind => atomic_kind_set(i)
2780 11189 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local)
2781 271498 : DO j = i, SIZE(atomic_kind_set)
2782 257688 : atomic_kind => atomic_kind_set(j)
2783 257688 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local)
2784 257688 : found = .FALSE.
2785 257688 : found_a = .FALSE.
2786 257688 : found_b = .FALSE.
2787 257688 : name_atm_a = name_atm_a_local
2788 257688 : name_atm_b = name_atm_b_local
2789 257688 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
2790 257688 : CALL uppercase(name_atm_a)
2791 257688 : CALL uppercase(name_atm_b)
2792 257688 : pot => potparm_nonbond14%pot(i, j)%pot
2793 :
2794 : ! loop over params from GROMOS
2795 257688 : IF (ASSOCIATED(gro_info%nonbond_a_14)) THEN
2796 540 : ii = 0
2797 540 : jj = 0
2798 1800 : DO k = 1, SIZE(gro_info%nonbond_a_14)
2799 1800 : IF (TRIM(name_atm_a) == TRIM(gro_info%nonbond_a_14(k))) THEN
2800 : ii = k
2801 : found_a = .TRUE.
2802 : EXIT
2803 : END IF
2804 : END DO
2805 2364 : DO k = 1, SIZE(gro_info%nonbond_a_14)
2806 2364 : IF (TRIM(name_atm_b) == TRIM(gro_info%nonbond_a_14(k))) THEN
2807 : jj = k
2808 : found_b = .TRUE.
2809 : EXIT
2810 : END IF
2811 : END DO
2812 540 : IF (ii /= 0 .AND. jj /= 0) THEN
2813 540 : CALL pair_potential_lj_create(pot%set(1)%lj)
2814 1080 : pot%type = lj_type
2815 540 : pot%at1 = name_atm_a
2816 540 : pot%at2 = name_atm_b
2817 540 : pot%set(1)%lj%epsilon = 1.0_dp
2818 540 : pot%set(1)%lj%sigma6 = gro_info%nonbond_c6_14(ii, jj)
2819 540 : pot%set(1)%lj%sigma12 = gro_info%nonbond_c12_14(ii, jj)
2820 540 : pot%rcutsq = (10.0_dp*bohr)**2
2821 540 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
2822 540 : found = .TRUE.
2823 : END IF
2824 : END IF
2825 :
2826 : ! Loop over params from CHARMM
2827 257688 : ii = 0
2828 257688 : jj = 0
2829 257688 : IF (ASSOCIATED(chm_info%nonbond_a_14)) THEN
2830 460416 : DO k = 1, SIZE(chm_info%nonbond_a_14)
2831 460416 : IF ((name_atm_a) == (chm_info%nonbond_a_14(k))) THEN
2832 11206 : ii = k
2833 11206 : rmin2_a = chm_info%nonbond_rmin2_14(k)
2834 11206 : epsilon_a = chm_info%nonbond_eps_14(k)
2835 11206 : found_a = .TRUE.
2836 : END IF
2837 : END DO
2838 460416 : DO k = 1, SIZE(chm_info%nonbond_a_14)
2839 460416 : IF ((name_atm_b) == (chm_info%nonbond_a_14(k))) THEN
2840 8888 : jj = k
2841 8888 : rmin2_b = chm_info%nonbond_rmin2_14(k)
2842 8888 : epsilon_b = chm_info%nonbond_eps_14(k)
2843 8888 : found_b = .TRUE.
2844 : END IF
2845 : END DO
2846 : END IF
2847 257688 : IF (ASSOCIATED(chm_info%nonbond_a)) THEN
2848 48329 : IF (.NOT. found_a) THEN
2849 1442245 : DO k = 1, SIZE(chm_info%nonbond_a)
2850 1442245 : IF ((name_atm_a) == (chm_info%nonbond_a(k))) THEN
2851 37057 : ii = k
2852 37057 : rmin2_a = chm_info%nonbond_rmin2(k)
2853 37057 : epsilon_a = chm_info%nonbond_eps(k)
2854 : END IF
2855 : END DO
2856 : END IF
2857 48329 : IF (.NOT. found_b) THEN
2858 1655155 : DO k = 1, SIZE(chm_info%nonbond_a)
2859 1655155 : IF ((name_atm_b) == (chm_info%nonbond_a(k))) THEN
2860 39423 : jj = k
2861 39423 : rmin2_b = chm_info%nonbond_rmin2(k)
2862 39423 : epsilon_b = chm_info%nonbond_eps(k)
2863 : END IF
2864 : END DO
2865 : END IF
2866 : END IF
2867 257688 : IF (ii /= 0 .AND. jj /= 0) THEN
2868 48263 : rmin = rmin2_a + rmin2_b
2869 : ! ABS to allow for mixing the two different sign conventions for epsilon
2870 48263 : epsilon0 = SQRT(ABS(epsilon_a*epsilon_b))
2871 48263 : CALL pair_potential_lj_create(pot%set(1)%lj)
2872 96526 : pot%type = lj_charmm_type
2873 48263 : pot%at1 = name_atm_a
2874 48263 : pot%at2 = name_atm_b
2875 48263 : pot%set(1)%lj%epsilon = epsilon0
2876 48263 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
2877 48263 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
2878 48263 : pot%rcutsq = (10.0_dp*bohr)**2
2879 48263 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
2880 48263 : found = .TRUE.
2881 : END IF
2882 :
2883 : ! Loop over params from AMBER
2884 257688 : IF (ASSOCIATED(amb_info%nonbond_a)) THEN
2885 199334 : ii = 0
2886 199334 : jj = 0
2887 199334 : IF (.NOT. found_a) THEN
2888 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
2889 45258092 : IF ((name_atm_a) == (amb_info%nonbond_a(k))) THEN
2890 199334 : ii = k
2891 199334 : rmin2_a = amb_info%nonbond_rmin2(k)
2892 199334 : epsilon_a = amb_info%nonbond_eps(k)
2893 : END IF
2894 : END DO
2895 : END IF
2896 199334 : IF (.NOT. found_b) THEN
2897 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
2898 45258092 : IF ((name_atm_b) == (amb_info%nonbond_a(k))) THEN
2899 199334 : jj = k
2900 199334 : rmin2_b = amb_info%nonbond_rmin2(k)
2901 199334 : epsilon_b = amb_info%nonbond_eps(k)
2902 : END IF
2903 : END DO
2904 : END IF
2905 199334 : IF (ii /= 0 .AND. jj /= 0) THEN
2906 199334 : rmin = rmin2_a + rmin2_b
2907 : ! ABS to allow for mixing the two different sign conventions for epsilon
2908 199334 : epsilon0 = SQRT(ABS(epsilon_a*epsilon_b))
2909 199334 : CALL pair_potential_lj_create(pot%set(1)%lj)
2910 398668 : pot%type = lj_charmm_type
2911 199334 : pot%at1 = name_atm_a
2912 199334 : pot%at2 = name_atm_b
2913 199334 : pot%set(1)%lj%epsilon = epsilon0
2914 199334 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
2915 199334 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
2916 199334 : pot%rcutsq = (10.0_dp*bohr)**2
2917 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, &
2918 199334 : name_atm_b)
2919 199334 : found = .TRUE.
2920 : END IF
2921 : END IF
2922 :
2923 : ! Always have the input param last to overwrite all the other ones
2924 257688 : IF (ASSOCIATED(inp_info%nonbonded14)) THEN
2925 12124 : DO k = 1, SIZE(inp_info%nonbonded14%pot)
2926 15265 : IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), &
2927 4817 : " with ", TRIM(inp_info%nonbonded14%pot(k)%pot%at1), &
2928 9634 : TRIM(inp_info%nonbonded14%pot(k)%pot%at2)
2929 : IF ((((name_atm_a) == (inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2930 10448 : ((name_atm_b) == (inp_info%nonbonded14%pot(k)%pot%at2))) .OR. &
2931 : (((name_atm_b) == (inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2932 1676 : ((name_atm_a) == (inp_info%nonbonded14%pot(k)%pot%at2)))) THEN
2933 1668 : IF (ff_type%multiple_potential) THEN
2934 0 : CALL pair_potential_single_add(inp_info%nonbonded14%pot(k)%pot, pot)
2935 0 : IF (found) THEN
2936 : CALL cp_warn(__LOCATION__, &
2937 : "Multiple ONFO declaration: "//TRIM(name_atm_a)// &
2938 0 : " and "//TRIM(name_atm_b)//" ADDING! ")
2939 : END IF
2940 0 : potparm_nonbond14%pot(i, j)%pot => pot
2941 0 : potparm_nonbond14%pot(j, i)%pot => pot
2942 : ELSE
2943 1668 : CALL pair_potential_single_copy(inp_info%nonbonded14%pot(k)%pot, pot)
2944 1668 : IF (found) THEN
2945 : CALL cp_warn(__LOCATION__, &
2946 : "Multiple ONFO declarations: "//TRIM(name_atm_a)// &
2947 0 : " and "//TRIM(name_atm_b)//" OVERWRITING! ")
2948 : END IF
2949 : END IF
2950 1668 : IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b)
2951 1668 : found = .TRUE.
2952 : END IF
2953 : END DO
2954 : END IF
2955 :
2956 : ! At the very end we offer the possibility to overwrite the parameters for QM/MM
2957 : ! nonbonded interactions
2958 257688 : IF (use_qmmm_ff) THEN
2959 252 : match_names = 0
2960 252 : IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1
2961 252 : IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1
2962 252 : IF (match_names == 1) THEN
2963 102 : IF (ASSOCIATED(qmmm_env%inp_info%nonbonded14)) THEN
2964 0 : DO k = 1, SIZE(qmmm_env%inp_info%nonbonded14%pot)
2965 : IF (debug_this_module) THEN
2966 : IF (iw > 0) THEN
2967 : WRITE (UNIT=iw, FMT="(T2,A)") &
2968 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
2969 : " with "//TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)//"-"// &
2970 : TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2)
2971 : END IF
2972 : END IF
2973 : IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2974 0 : ((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2))) .OR. &
2975 : (((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. &
2976 0 : ((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2)))) THEN
2977 0 : IF (qmmm_env%multiple_potential) THEN
2978 0 : CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded14%pot(k)%pot, pot)
2979 0 : IF (found) THEN
2980 : CALL cp_warn(__LOCATION__, &
2981 : "Multiple ONFO declaration: "//TRIM(name_atm_a)// &
2982 0 : " and "//TRIM(name_atm_b)//" Adding QM/MM forcefield specifications")
2983 : END IF
2984 0 : potparm_nonbond14%pot(i, j)%pot => pot
2985 0 : potparm_nonbond14%pot(j, i)%pot => pot
2986 : ELSE
2987 0 : CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded14%pot(k)%pot, pot)
2988 0 : IF (found) THEN
2989 : CALL cp_warn(__LOCATION__, &
2990 : "Multiple ONFO declaration: "//TRIM(name_atm_a)// &
2991 0 : " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! ")
2992 : END IF
2993 : END IF
2994 0 : IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), &
2995 0 : " ", TRIM(name_atm_b)
2996 0 : found = .TRUE.
2997 : END IF
2998 : END DO
2999 : END IF
3000 : END IF
3001 : END IF
3002 :
3003 257688 : IF (.NOT. found) THEN
3004 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
3005 : atm2=TRIM(name_atm_b), &
3006 : type_name="Spline_Bond_Env", &
3007 7883 : array=Ainfo)
3008 7883 : CALL pair_potential_single_clean(pot)
3009 15766 : pot%type = nn_type
3010 7883 : pot%at1 = name_atm_a
3011 7883 : pot%at2 = name_atm_b
3012 : END IF
3013 :
3014 : ! If defined global RCUT let's use it
3015 257688 : IF (ff_type%rcut_nb > 0.0_dp) THEN
3016 26948 : pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb
3017 : END IF
3018 :
3019 : ! Cutoff is defined always as the maximum between the FF and Ewald
3020 257688 : pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut)
3021 268877 : IF (only_qm) THEN
3022 11786 : CALL pair_potential_single_clean(pot)
3023 : END IF
3024 :
3025 : END DO ! atom kind j
3026 :
3027 : END DO ! atom kind i
3028 :
3029 2621 : CALL timestop(handle2)
3030 :
3031 2621 : END SUBROUTINE force_field_pack_nonbond14
3032 :
3033 : ! **************************************************************************************************
3034 : !> \brief Assign input and potential info to potparm_nonbond
3035 : !> \param atomic_kind_set ...
3036 : !> \param ff_type ...
3037 : !> \param qmmm_env ...
3038 : !> \param fatal ...
3039 : !> \param iw ...
3040 : !> \param Ainfo ...
3041 : !> \param chm_info ...
3042 : !> \param inp_info ...
3043 : !> \param gro_info ...
3044 : !> \param amb_info ...
3045 : !> \param potparm_nonbond ...
3046 : !> \param ewald_env ...
3047 : ! **************************************************************************************************
3048 2621 : SUBROUTINE force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, fatal, &
3049 : iw, Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond, &
3050 : ewald_env)
3051 :
3052 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3053 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
3054 : TYPE(qmmm_env_mm_type), POINTER :: qmmm_env
3055 : LOGICAL :: fatal
3056 : INTEGER :: iw
3057 : CHARACTER(LEN=default_string_length), &
3058 : DIMENSION(:), POINTER :: Ainfo
3059 : TYPE(charmm_info_type), POINTER :: chm_info
3060 : TYPE(input_info_type), POINTER :: inp_info
3061 : TYPE(gromos_info_type), POINTER :: gro_info
3062 : TYPE(amber_info_type), POINTER :: amb_info
3063 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond
3064 : TYPE(ewald_environment_type), POINTER :: ewald_env
3065 :
3066 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond'
3067 :
3068 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, &
3069 : name_atm_b, name_atm_b_local
3070 : INTEGER :: handle2, i, ii, j, jj, k, match_names
3071 : LOGICAL :: found, is_a_shell, is_b_shell, only_qm, &
3072 : use_qmmm_ff
3073 : REAL(KIND=dp) :: epsilon0, ewald_rcut, rmin
3074 : TYPE(atomic_kind_type), POINTER :: atomic_kind
3075 : TYPE(pair_potential_single_type), POINTER :: pot
3076 :
3077 2621 : CALL timeset(routineN, handle2)
3078 :
3079 2621 : use_qmmm_ff = qmmm_env%use_qmmm_ff
3080 2621 : NULLIFY (pot)
3081 :
3082 2621 : IF (iw > 0) THEN
3083 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
3084 241 : "FORCEFIELD| Checking for nonbonded terms"
3085 : END IF
3086 :
3087 2621 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
3088 2621 : CALL pair_potential_pp_create(potparm_nonbond, SIZE(atomic_kind_set))
3089 :
3090 13810 : DO i = 1, SIZE(atomic_kind_set)
3091 :
3092 11189 : atomic_kind => atomic_kind_set(i)
3093 :
3094 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local, &
3095 11189 : shell_active=is_a_shell)
3096 :
3097 271498 : DO j = i, SIZE(atomic_kind_set)
3098 :
3099 257688 : atomic_kind => atomic_kind_set(j)
3100 :
3101 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local, &
3102 257688 : shell_active=is_b_shell)
3103 :
3104 257688 : found = .FALSE.
3105 :
3106 257688 : name_atm_a = name_atm_a_local
3107 257688 : name_atm_b = name_atm_b_local
3108 257688 : only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
3109 257688 : CALL uppercase(name_atm_a)
3110 257688 : CALL uppercase(name_atm_b)
3111 257688 : pot => potparm_nonbond%pot(i, j)%pot
3112 :
3113 257688 : IF (iw > 0) THEN
3114 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
3115 : "FORCEFIELD| Checking for nonbonded terms between the atomic kinds "// &
3116 6444 : TRIM(name_atm_a)//" and "//TRIM(name_atm_b)
3117 : END IF
3118 :
3119 : ! Loop over params from GROMOS
3120 257688 : IF (ASSOCIATED(gro_info%nonbond_a)) THEN
3121 540 : ii = 0
3122 540 : jj = 0
3123 1800 : DO k = 1, SIZE(gro_info%nonbond_a)
3124 1800 : IF (TRIM(name_atm_a) == TRIM(gro_info%nonbond_a(k))) THEN
3125 : ii = k
3126 : EXIT
3127 : END IF
3128 : END DO
3129 2364 : DO k = 1, SIZE(gro_info%nonbond_a)
3130 2364 : IF (TRIM(name_atm_b) == TRIM(gro_info%nonbond_a(k))) THEN
3131 : jj = k
3132 : EXIT
3133 : END IF
3134 : END DO
3135 :
3136 540 : IF (ii /= 0 .AND. jj /= 0) THEN
3137 540 : CALL pair_potential_lj_create(pot%set(1)%lj)
3138 1080 : pot%type = lj_type
3139 540 : pot%at1 = name_atm_a
3140 540 : pot%at2 = name_atm_b
3141 540 : pot%set(1)%lj%epsilon = 1.0_dp
3142 540 : pot%set(1)%lj%sigma6 = gro_info%nonbond_c6(ii, jj)
3143 540 : pot%set(1)%lj%sigma12 = gro_info%nonbond_c12(ii, jj)
3144 540 : pot%rcutsq = (10.0_dp*bohr)**2
3145 540 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
3146 540 : found = .TRUE.
3147 : END IF
3148 : END IF
3149 :
3150 : ! Loop over params from CHARMM
3151 257688 : IF (ASSOCIATED(chm_info%nonbond_a)) THEN
3152 48329 : ii = 0
3153 48329 : jj = 0
3154 2286557 : DO k = 1, SIZE(chm_info%nonbond_a)
3155 2286557 : IF ((name_atm_a) == (chm_info%nonbond_a(k))) THEN
3156 48263 : ii = k
3157 : END IF
3158 : END DO
3159 2286557 : DO k = 1, SIZE(chm_info%nonbond_a)
3160 2286557 : IF ((name_atm_b) == (chm_info%nonbond_a(k))) THEN
3161 48311 : jj = k
3162 : END IF
3163 : END DO
3164 :
3165 48329 : IF (ii /= 0 .AND. jj /= 0) THEN
3166 48263 : rmin = chm_info%nonbond_rmin2(ii) + chm_info%nonbond_rmin2(jj)
3167 : epsilon0 = SQRT(chm_info%nonbond_eps(ii)* &
3168 48263 : chm_info%nonbond_eps(jj))
3169 48263 : CALL pair_potential_lj_create(pot%set(1)%lj)
3170 96526 : pot%type = lj_charmm_type
3171 48263 : pot%at1 = name_atm_a
3172 48263 : pot%at2 = name_atm_b
3173 48263 : pot%set(1)%lj%epsilon = epsilon0
3174 48263 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
3175 48263 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
3176 48263 : pot%rcutsq = (10.0_dp*bohr)**2
3177 48263 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
3178 48263 : found = .TRUE.
3179 : END IF
3180 : END IF
3181 :
3182 : ! Loop over params from AMBER
3183 257688 : IF (ASSOCIATED(amb_info%nonbond_a)) THEN
3184 199334 : ii = 0
3185 199334 : jj = 0
3186 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
3187 45258092 : IF ((name_atm_a) == (amb_info%nonbond_a(k))) THEN
3188 199334 : ii = k
3189 : END IF
3190 : END DO
3191 45258092 : DO k = 1, SIZE(amb_info%nonbond_a)
3192 45258092 : IF ((name_atm_b) == (amb_info%nonbond_a(k))) THEN
3193 199334 : jj = k
3194 : END IF
3195 : END DO
3196 :
3197 199334 : IF (ii /= 0 .AND. jj /= 0) THEN
3198 199334 : rmin = amb_info%nonbond_rmin2(ii) + amb_info%nonbond_rmin2(jj)
3199 199334 : epsilon0 = SQRT(amb_info%nonbond_eps(ii)*amb_info%nonbond_eps(jj))
3200 199334 : CALL pair_potential_lj_create(pot%set(1)%lj)
3201 398668 : pot%type = lj_charmm_type
3202 199334 : pot%at1 = name_atm_a
3203 199334 : pot%at2 = name_atm_b
3204 199334 : pot%set(1)%lj%epsilon = epsilon0
3205 199334 : pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6
3206 199334 : pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12
3207 199334 : pot%rcutsq = (10.0_dp*bohr)**2
3208 199334 : CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b)
3209 199334 : found = .TRUE.
3210 : END IF
3211 : END IF
3212 :
3213 : ! Always have the input param last to overwrite all the other ones
3214 257688 : IF (ASSOCIATED(inp_info%nonbonded)) THEN
3215 53070 : DO k = 1, SIZE(inp_info%nonbonded%pot)
3216 43439 : IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .OR. &
3217 : (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE
3218 : IF (debug_this_module) THEN
3219 : IF (iw > 0) THEN
3220 : WRITE (UNIT=iw, FMT="(T2,A)") &
3221 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3222 : " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// &
3223 : TRIM(inp_info%nonbonded%pot(k)%pot%at2)
3224 : END IF
3225 : END IF
3226 : IF ((((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3227 43437 : ((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at2))) .OR. &
3228 : (((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3229 9631 : ((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at2)))) THEN
3230 9410 : IF (iw > 0) THEN
3231 : WRITE (UNIT=iw, FMT="(T2,A)") &
3232 : "FORCEFIELD| Found nonbonded term "// &
3233 935 : TRIM(name_atm_a)//"-"//TRIM(name_atm_b)
3234 : END IF
3235 9410 : IF (ff_type%multiple_potential) THEN
3236 38 : CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot)
3237 38 : IF (found) THEN
3238 : CALL cp_warn(__LOCATION__, &
3239 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3240 8 : "-"//TRIM(name_atm_b)//" -> ADDING")
3241 : END IF
3242 38 : potparm_nonbond%pot(i, j)%pot => pot
3243 38 : potparm_nonbond%pot(j, i)%pot => pot
3244 : ELSE
3245 9372 : CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot)
3246 9372 : IF (found) THEN
3247 : CALL cp_warn(__LOCATION__, &
3248 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3249 8 : "-"//TRIM(name_atm_b)//" -> OVERWRITING")
3250 : END IF
3251 : END IF
3252 9410 : found = .TRUE.
3253 : END IF
3254 : END DO
3255 :
3256 : ! Check for wildcards for one of the two types (if not associated yet)
3257 9631 : IF (.NOT. found) THEN
3258 590 : DO k = 1, SIZE(inp_info%nonbonded%pot)
3259 433 : IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .EQV. &
3260 : (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE
3261 : IF (debug_this_module) THEN
3262 : IF (iw > 0) THEN
3263 : WRITE (UNIT=iw, FMT="(T2,A)") &
3264 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3265 : " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// &
3266 : TRIM(inp_info%nonbonded%pot(k)%pot%at2)
3267 : END IF
3268 : END IF
3269 : IF ((name_atm_a == inp_info%nonbonded%pot(k)%pot%at1) .OR. &
3270 : (name_atm_b == inp_info%nonbonded%pot(k)%pot%at2) .OR. &
3271 0 : (name_atm_b == inp_info%nonbonded%pot(k)%pot%at1) .OR. &
3272 157 : (name_atm_a == inp_info%nonbonded%pot(k)%pot%at2)) THEN
3273 0 : IF (iw > 0) THEN
3274 : WRITE (UNIT=iw, FMT="(T2,A)") &
3275 : "FORCEFIELD| Found one wildcard for "// &
3276 0 : TRIM(name_atm_a)//"-"//TRIM(name_atm_b)
3277 : END IF
3278 0 : IF (ff_type%multiple_potential) THEN
3279 0 : CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot)
3280 0 : IF (found) THEN
3281 : CALL cp_warn(__LOCATION__, &
3282 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3283 0 : "-"//TRIM(name_atm_b)//" -> ADDING")
3284 : END IF
3285 0 : potparm_nonbond%pot(i, j)%pot => pot
3286 0 : potparm_nonbond%pot(j, i)%pot => pot
3287 : ELSE
3288 0 : CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot)
3289 0 : IF (found) THEN
3290 : CALL cp_warn(__LOCATION__, &
3291 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3292 0 : "-"//TRIM(name_atm_b)//" -> OVERWRITING")
3293 : END IF
3294 : END IF
3295 0 : found = .TRUE.
3296 : END IF
3297 : END DO
3298 : END IF
3299 :
3300 : ! Check for wildcards for both types (if not associated yet)
3301 9631 : IF (.NOT. found) THEN
3302 590 : DO k = 1, SIZE(inp_info%nonbonded%pot)
3303 433 : IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) /= "*") .OR. &
3304 : (TRIM(inp_info%nonbonded%pot(k)%pot%at2) /= "*")) CYCLE
3305 : IF (debug_this_module) THEN
3306 : IF (iw > 0) THEN
3307 : WRITE (UNIT=iw, FMT="(T2,A)") &
3308 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3309 : " with "//TRIM(inp_info%nonbonded%pot(k)%pot%at1)//"-"// &
3310 : TRIM(inp_info%nonbonded%pot(k)%pot%at2)
3311 : END IF
3312 : END IF
3313 2 : IF (iw > 0) THEN
3314 : WRITE (UNIT=iw, FMT="(T2,A)") &
3315 : "FORCEFIELD| Found wildcards for both "// &
3316 0 : TRIM(name_atm_a)//" and "//TRIM(name_atm_b)
3317 : END IF
3318 2 : IF (ff_type%multiple_potential) THEN
3319 0 : CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot)
3320 0 : IF (found) THEN
3321 : CALL cp_warn(__LOCATION__, &
3322 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3323 0 : " - "//TRIM(name_atm_b)//" -> ADDING")
3324 : END IF
3325 0 : potparm_nonbond%pot(i, j)%pot => pot
3326 0 : potparm_nonbond%pot(j, i)%pot => pot
3327 : ELSE
3328 2 : CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot)
3329 2 : IF (found) THEN
3330 : CALL cp_warn(__LOCATION__, &
3331 : "Multiple NONBONDED declarations "//TRIM(name_atm_a)// &
3332 0 : " - "//TRIM(name_atm_b)//" -> OVERWRITING")
3333 : END IF
3334 : END IF
3335 590 : found = .TRUE.
3336 : END DO
3337 : END IF
3338 : END IF
3339 :
3340 : ! At the very end we offer the possibility to overwrite the parameters for QM/MM
3341 : ! nonbonded interactions
3342 257688 : IF (use_qmmm_ff) THEN
3343 252 : match_names = 0
3344 252 : IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1
3345 252 : IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1
3346 252 : IF (match_names == 1) THEN
3347 102 : IF (ASSOCIATED(qmmm_env%inp_info%nonbonded)) THEN
3348 276 : DO k = 1, SIZE(qmmm_env%inp_info%nonbonded%pot)
3349 : IF (debug_this_module) THEN
3350 : IF (iw > 0) THEN
3351 : WRITE (UNIT=iw, FMT="(T2,A)") &
3352 : "FORCEFIELD| Testing "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)// &
3353 : " with "//TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1), &
3354 : TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2)
3355 : END IF
3356 : END IF
3357 : IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3358 174 : ((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2))) .OR. &
3359 : (((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. &
3360 102 : ((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2)))) THEN
3361 20 : IF (iw > 0) THEN
3362 : WRITE (UNIT=iw, FMT="(T2,A)") &
3363 11 : "FORCEFIELD| Found "//TRIM(name_atm_a)//"-"//TRIM(name_atm_b)//" (QM/MM)"
3364 : END IF
3365 20 : IF (qmmm_env%multiple_potential) THEN
3366 0 : CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot)
3367 0 : IF (found) THEN
3368 : CALL cp_warn(__LOCATION__, &
3369 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3370 0 : " and "//TRIM(name_atm_b)//" -> ADDING QM/MM forcefield specifications")
3371 : END IF
3372 0 : potparm_nonbond%pot(i, j)%pot => pot
3373 0 : potparm_nonbond%pot(j, i)%pot => pot
3374 : ELSE
3375 20 : CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot)
3376 20 : IF (found) THEN
3377 : CALL cp_warn(__LOCATION__, &
3378 : "Multiple NONBONDED declarations for "//TRIM(name_atm_a)// &
3379 2 : " and "//TRIM(name_atm_b)//" -> OVERWRITING QM/MM forcefield specifications")
3380 : END IF
3381 : END IF
3382 20 : found = .TRUE.
3383 : END IF
3384 : END DO
3385 : END IF
3386 : END IF
3387 : END IF
3388 :
3389 257688 : IF (.NOT. found) THEN
3390 : CALL store_FF_missing_par(atm1=TRIM(name_atm_a), &
3391 : atm2=TRIM(name_atm_b), &
3392 : type_name="Spline_Non_Bond_Env", &
3393 : fatal=fatal, &
3394 137 : array=Ainfo)
3395 : END IF
3396 :
3397 : ! If defined global RCUT let's use it
3398 257688 : IF (ff_type%rcut_nb > 0.0_dp) THEN
3399 26948 : pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb
3400 : END IF
3401 :
3402 : ! Cutoff is defined always as the maximum between the FF and Ewald
3403 257688 : pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut)
3404 : ! Set the shell type
3405 257688 : IF ((is_a_shell .AND. .NOT. is_b_shell) .OR. (is_b_shell .AND. .NOT. is_a_shell)) THEN
3406 56 : pot%shell_type = nosh_sh
3407 257632 : ELSE IF (is_a_shell .AND. is_b_shell) THEN
3408 624 : pot%shell_type = sh_sh
3409 : ELSE
3410 257008 : pot%shell_type = nosh_nosh
3411 : END IF
3412 :
3413 526565 : IF (only_qm) THEN
3414 11786 : CALL pair_potential_single_clean(pot)
3415 : END IF
3416 :
3417 : END DO ! jkind
3418 :
3419 : END DO ! ikind
3420 :
3421 2621 : CALL timestop(handle2)
3422 :
3423 2621 : END SUBROUTINE force_field_pack_nonbond
3424 :
3425 : ! **************************************************************************************************
3426 : !> \brief create the pair potential spline environment
3427 : !> \param atomic_kind_set ...
3428 : !> \param ff_type ...
3429 : !> \param iw2 ...
3430 : !> \param iw3 ...
3431 : !> \param iw4 ...
3432 : !> \param potparm ...
3433 : !> \param do_zbl ...
3434 : !> \param nonbonded_type ...
3435 : ! **************************************************************************************************
3436 5242 : SUBROUTINE force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, &
3437 : potparm, do_zbl, nonbonded_type)
3438 :
3439 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3440 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
3441 : INTEGER :: iw2, iw3, iw4
3442 : TYPE(pair_potential_pp_type), POINTER :: potparm
3443 : LOGICAL, INTENT(IN) :: do_zbl
3444 : CHARACTER(LEN=*), INTENT(IN) :: nonbonded_type
3445 :
3446 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_splines'
3447 :
3448 : INTEGER :: handle2, ikind, jkind, n
3449 5242 : TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spl_p
3450 : TYPE(spline_environment_type), POINTER :: spline_env
3451 :
3452 5242 : CALL timeset(routineN, handle2)
3453 :
3454 5242 : IF (iw2 > 0) THEN
3455 : WRITE (UNIT=iw2, FMT="(/,T2,A)") &
3456 482 : "FORCEFIELD| Splining nonbonded terms"
3457 : END IF
3458 :
3459 : ! Figure out which nonbonded interactions happen to be identical, and
3460 : ! prepare storage for these, avoiding duplicates.
3461 5242 : NULLIFY (spline_env)
3462 : CALL get_nonbond_storage(spline_env, potparm, atomic_kind_set, &
3463 5242 : do_zbl, shift_cutoff=ff_type%shift_cutoff)
3464 : ! Effectively compute the spline data
3465 : CALL spline_nonbond_control(spline_env, potparm, &
3466 : atomic_kind_set, eps_spline=ff_type%eps_spline, &
3467 : max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb, &
3468 : emax_spline=ff_type%emax_spline, npoints=ff_type%npoints, &
3469 : iw=iw2, iw2=iw3, iw3=iw4, &
3470 : do_zbl=do_zbl, shift_cutoff=ff_type%shift_cutoff, &
3471 5242 : nonbonded_type=nonbonded_type)
3472 : ! Let the pointers on potparm point to the splines generated in
3473 : ! spline_nonbond_control
3474 27620 : DO ikind = 1, SIZE(potparm%pot, 1)
3475 542996 : DO jkind = ikind, SIZE(potparm%pot, 2)
3476 515376 : n = spline_env%spltab(ikind, jkind)
3477 515376 : spl_p => spline_env%spl_pp(n)%spl_p
3478 515376 : CALL spline_data_p_retain(spl_p)
3479 515376 : CALL spline_data_p_release(potparm%pot(ikind, jkind)%pot%pair_spline_data)
3480 537754 : potparm%pot(ikind, jkind)%pot%pair_spline_data => spl_p
3481 : END DO
3482 : END DO
3483 5242 : CALL spline_env_release(spline_env)
3484 5242 : DEALLOCATE (spline_env)
3485 : NULLIFY (spline_env)
3486 :
3487 5242 : IF (iw2 > 0) THEN
3488 : WRITE (UNIT=iw2, FMT="(/,T2,A)") &
3489 482 : "FORCEFIELD| Splining done"
3490 : END IF
3491 :
3492 5242 : CALL timestop(handle2)
3493 :
3494 5242 : END SUBROUTINE force_field_pack_splines
3495 :
3496 : ! **************************************************************************************************
3497 : !> \brief Compute the electrostatic interaction cutoffs
3498 : !> \param atomic_kind_set ...
3499 : !> \param ff_type ...
3500 : !> \param potparm_nonbond ...
3501 : !> \param ewald_env ...
3502 : !> \param iw ...
3503 : !> \author Toon.Verstraelen@gmail.com
3504 : ! **************************************************************************************************
3505 2637 : SUBROUTINE force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, ewald_env, iw)
3506 :
3507 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3508 : TYPE(force_field_type), INTENT(IN) :: ff_type
3509 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond
3510 : TYPE(ewald_environment_type), POINTER :: ewald_env
3511 : INTEGER, INTENT(IN) :: iw
3512 :
3513 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_eicut'
3514 :
3515 : INTEGER :: ewald_type, handle, i1, i2, nkinds
3516 : REAL(KIND=dp) :: alpha, beta, mm_radius1, mm_radius2, &
3517 : rcut2, rcut2_ewald, tmp
3518 2637 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs
3519 : TYPE(atomic_kind_type), POINTER :: atomic_kind
3520 :
3521 2637 : CALL timeset(routineN, handle)
3522 :
3523 2637 : IF (iw > 0) THEN
3524 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
3525 242 : "FORCEFIELD| Computing the electrostatic interactions cutoffs"
3526 : END IF
3527 :
3528 2637 : tmp = 0.0_dp
3529 2637 : nkinds = SIZE(atomic_kind_set)
3530 :
3531 : ! Allocate the array with interaction cutoffs for the electrostatics, used
3532 : ! to make the electrostatic interaction continuous at ewald_env%rcut
3533 10548 : ALLOCATE (interaction_cutoffs(3, nkinds, nkinds))
3534 2032768 : interaction_cutoffs = 0.0_dp
3535 :
3536 : ! Compute the interaction cutoff if SHIFT_CUTOFF is active
3537 2637 : IF (ff_type%shift_cutoff) THEN
3538 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type, &
3539 2485 : rcut=rcut2_ewald)
3540 2485 : rcut2_ewald = rcut2_ewald*rcut2_ewald
3541 11708 : DO i1 = 1, nkinds
3542 9223 : atomic_kind => atomic_kind_set(i1)
3543 9223 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius1)
3544 118045 : DO i2 = 1, nkinds
3545 106337 : rcut2 = rcut2_ewald
3546 106337 : IF (ASSOCIATED(potparm_nonbond)) THEN
3547 105807 : rcut2 = MAX(potparm_nonbond%pot(i1, i2)%pot%rcutsq, rcut2_ewald)
3548 : END IF
3549 115560 : IF (rcut2 > 0) THEN
3550 103025 : atomic_kind => atomic_kind_set(i2)
3551 103025 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius2)
3552 : ! cutoff for core-core
3553 : interaction_cutoffs(1, i1, i2) = potential_coulomb(rcut2, tmp, &
3554 103025 : 1.0_dp, ewald_type, alpha, 0.0_dp, 0.0_dp)
3555 : ! cutoff for core-shell, core-ion, shell-core or ion-core
3556 103025 : IF (mm_radius1 > 0.0_dp) THEN
3557 676 : beta = sqrthalf/mm_radius1
3558 : ELSE
3559 102349 : beta = 0.0_dp
3560 : END IF
3561 : interaction_cutoffs(2, i1, i2) = potential_coulomb(rcut2, tmp, &
3562 103025 : 1.0_dp, ewald_type, alpha, beta, 0.0_dp)
3563 : ! cutoff for shell-shell or ion-ion
3564 103025 : IF (mm_radius1 + mm_radius2 > 0.0_dp) THEN
3565 698 : beta = sqrthalf/SQRT(mm_radius1*mm_radius1 + mm_radius2*mm_radius2)
3566 : ELSE
3567 102327 : beta = 0.0_dp
3568 : END IF
3569 : interaction_cutoffs(3, i1, i2) = potential_coulomb(rcut2, tmp, &
3570 103025 : 1.0_dp, ewald_type, alpha, beta, 0.0_dp)
3571 : END IF
3572 : END DO
3573 : END DO
3574 : END IF
3575 :
3576 2637 : CALL ewald_env_set(ewald_env, interaction_cutoffs=interaction_cutoffs)
3577 :
3578 2637 : CALL timestop(handle)
3579 :
3580 2637 : END SUBROUTINE force_field_pack_eicut
3581 :
3582 : ! **************************************************************************************************
3583 : !> \brief Issues on screen a warning when repetitions are present in the
3584 : !> definition of the forcefield
3585 : !> \param found ...
3586 : !> \param tag_label ...
3587 : !> \param name_atm_a ...
3588 : !> \param name_atm_b ...
3589 : !> \param name_atm_c ...
3590 : !> \param name_atm_d ...
3591 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
3592 : ! **************************************************************************************************
3593 783537 : SUBROUTINE issue_duplications(found, tag_label, name_atm_a, name_atm_b, &
3594 : name_atm_c, name_atm_d)
3595 :
3596 : LOGICAL, INTENT(IN) :: found
3597 : CHARACTER(LEN=*), INTENT(IN) :: tag_label, name_atm_a
3598 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name_atm_b, name_atm_c, name_atm_d
3599 :
3600 : CHARACTER(LEN=default_string_length) :: item
3601 :
3602 783537 : item = "("//TRIM(name_atm_a)
3603 783537 : IF (PRESENT(name_atm_b)) THEN
3604 775701 : item = TRIM(item)//", "//TRIM(name_atm_b)
3605 : END IF
3606 783537 : IF (PRESENT(name_atm_c)) THEN
3607 164596 : item = TRIM(item)//", "//TRIM(name_atm_c)
3608 : END IF
3609 783537 : IF (PRESENT(name_atm_d)) THEN
3610 3418 : item = TRIM(item)//", "//TRIM(name_atm_d)
3611 : END IF
3612 783537 : item = TRIM(item)//")"
3613 783537 : IF (found) THEN
3614 1678 : CPWARN("Found multiple "//TRIM(tag_label)//" terms for "//TRIM(item)//" -> OVERWRITING")
3615 : END IF
3616 :
3617 783537 : END SUBROUTINE issue_duplications
3618 :
3619 : ! **************************************************************************************************
3620 : !> \brief Store informations on possible missing ForceFields parameters
3621 : !> \param atm1 ...
3622 : !> \param atm2 ...
3623 : !> \param atm3 ...
3624 : !> \param atm4 ...
3625 : !> \param type_name ...
3626 : !> \param fatal ...
3627 : !> \param array ...
3628 : ! **************************************************************************************************
3629 171120 : SUBROUTINE store_FF_missing_par(atm1, atm2, atm3, atm4, type_name, fatal, array)
3630 : CHARACTER(LEN=*), INTENT(IN) :: atm1
3631 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: atm2, atm3, atm4
3632 : CHARACTER(LEN=*), INTENT(IN) :: type_name
3633 : LOGICAL, INTENT(INOUT), OPTIONAL :: fatal
3634 : CHARACTER(LEN=default_string_length), &
3635 : DIMENSION(:), POINTER :: array
3636 :
3637 : CHARACTER(LEN=10) :: sfmt
3638 : CHARACTER(LEN=9) :: my_atm1, my_atm2, my_atm3, my_atm4
3639 : CHARACTER(LEN=default_path_length) :: my_format
3640 : INTEGER :: fmt, i, nsize
3641 : LOGICAL :: found
3642 :
3643 171120 : nsize = 0
3644 171120 : fmt = 1
3645 : my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)// &
3646 171120 : '",T40,"(",A9,")")'
3647 171120 : IF (PRESENT(atm2)) fmt = fmt + 1
3648 171120 : IF (PRESENT(atm3)) fmt = fmt + 1
3649 171120 : IF (PRESENT(atm4)) fmt = fmt + 1
3650 171120 : CALL integer_to_string(fmt - 1, sfmt)
3651 171120 : IF (fmt > 1) THEN
3652 : my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)// &
3653 171108 : '",T40,"(",A9,'//TRIM(sfmt)//'(",",A9),")")'
3654 : END IF
3655 171120 : IF (PRESENT(fatal)) fatal = .TRUE.
3656 : ! Check for previous already stored equal force fields
3657 171120 : IF (ASSOCIATED(array)) nsize = SIZE(array)
3658 171120 : found = .FALSE.
3659 171120 : IF (nsize >= 1) THEN
3660 19488308 : DO i = 1, nsize
3661 8 : SELECT CASE (type_name)
3662 : CASE ("Bond")
3663 8 : IF (INDEX(array(i) (21:39), "Bond") == 0) CYCLE
3664 8 : my_atm1 = array(i) (41:49)
3665 8 : my_atm2 = array(i) (51:59)
3666 8 : CALL compress(my_atm1, .TRUE.)
3667 8 : CALL compress(my_atm2, .TRUE.)
3668 8 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. &
3669 8 : ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .TRUE.
3670 : CASE ("Angle")
3671 8 : IF (INDEX(array(i) (21:39), "Angle") == 0) CYCLE
3672 0 : my_atm1 = array(i) (41:49)
3673 0 : my_atm2 = array(i) (51:59)
3674 0 : my_atm3 = array(i) (61:69)
3675 0 : CALL compress(my_atm1, .TRUE.)
3676 0 : CALL compress(my_atm2, .TRUE.)
3677 0 : CALL compress(my_atm3, .TRUE.)
3678 0 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. &
3679 18206414 : ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) THEN
3680 : found = .TRUE.
3681 : END IF
3682 : CASE ("Urey-Bradley")
3683 18206414 : IF (INDEX(array(i) (21:39), "Urey-Bradley") == 0) CYCLE
3684 18206414 : my_atm1 = array(i) (41:49)
3685 18206414 : my_atm2 = array(i) (51:59)
3686 18206414 : my_atm3 = array(i) (61:69)
3687 18206414 : CALL compress(my_atm1, .TRUE.)
3688 18206414 : CALL compress(my_atm2, .TRUE.)
3689 18206414 : CALL compress(my_atm3, .TRUE.)
3690 18206414 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. &
3691 607150 : ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) THEN
3692 : found = .TRUE.
3693 : END IF
3694 : CASE ("Torsion")
3695 607150 : IF (INDEX(array(i) (21:39), "Torsion") == 0) CYCLE
3696 198764 : my_atm1 = array(i) (41:49)
3697 198764 : my_atm2 = array(i) (51:59)
3698 198764 : my_atm3 = array(i) (61:69)
3699 198764 : my_atm4 = array(i) (71:79)
3700 198764 : CALL compress(my_atm1, .TRUE.)
3701 198764 : CALL compress(my_atm2, .TRUE.)
3702 198764 : CALL compress(my_atm3, .TRUE.)
3703 198764 : CALL compress(my_atm4, .TRUE.)
3704 198764 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3705 154212 : ((atm1 == my_atm4) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm1))) THEN
3706 : found = .TRUE.
3707 : END IF
3708 : CASE ("Improper")
3709 154212 : IF (INDEX(array(i) (21:39), "Improper") == 0) CYCLE
3710 9684 : my_atm1 = array(i) (41:49)
3711 9684 : my_atm2 = array(i) (51:59)
3712 9684 : my_atm3 = array(i) (61:69)
3713 9684 : my_atm4 = array(i) (71:79)
3714 9684 : CALL compress(my_atm1, .TRUE.)
3715 9684 : CALL compress(my_atm2, .TRUE.)
3716 9684 : CALL compress(my_atm3, .TRUE.)
3717 9684 : CALL compress(my_atm4, .TRUE.)
3718 : 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)) .OR. &
3720 : ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3)) .OR. &
3721 : ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm2)) .OR. &
3722 9684 : ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm3)) .OR. &
3723 483920 : ((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3))) THEN
3724 : found = .TRUE.
3725 : END IF
3726 :
3727 : CASE ("Out of plane bend")
3728 483920 : IF (INDEX(array(i) (21:39), "Out of plane bend") == 0) CYCLE
3729 27416 : my_atm1 = array(i) (41:49)
3730 27416 : my_atm2 = array(i) (51:59)
3731 27416 : my_atm3 = array(i) (61:69)
3732 27416 : my_atm4 = array(i) (71:79)
3733 27416 : CALL compress(my_atm1, .TRUE.)
3734 27416 : CALL compress(my_atm2, .TRUE.)
3735 27416 : CALL compress(my_atm3, .TRUE.)
3736 27416 : CALL compress(my_atm4, .TRUE.)
3737 27416 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3738 8 : ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4))) THEN
3739 : found = .TRUE.
3740 : END IF
3741 :
3742 : CASE ("Charge")
3743 8 : IF (INDEX(array(i) (21:39), "Charge") == 0) CYCLE
3744 8 : my_atm1 = array(i) (41:49)
3745 8 : CALL compress(my_atm1, .TRUE.)
3746 8 : IF (atm1 == my_atm1) found = .TRUE.
3747 : CASE ("Spline_Bond_Env", "Spline_Non_Bond_Env")
3748 18778 : IF (INDEX(array(i) (21:39), "Spline_") == 0) CYCLE
3749 6539 : fmt = 0
3750 6539 : my_atm1 = array(i) (41:49)
3751 6539 : my_atm2 = array(i) (51:59)
3752 6539 : CALL compress(my_atm1, .TRUE.)
3753 6539 : CALL compress(my_atm2, .TRUE.)
3754 6539 : IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. &
3755 0 : ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .TRUE.
3756 : CASE DEFAULT
3757 : ! Should never reach this point
3758 : CALL cp_abort(__LOCATION__, &
3759 : "<Bond>, <Angle>, <Urey-Bradley>, <Torsion>, "// &
3760 : "<Improper>, <Out of plane bend>, <Charge>, "// &
3761 : "<Spline_Bond_Env>, <Spline_Non_Bond_Env> are "// &
3762 : "supported as the <type_name> for store_FF_missing_par, "// &
3763 : "found unknown option "// &
3764 19470498 : "<"//TRIM(type_name)//">")
3765 : END SELECT
3766 17810 : IF (found) EXIT
3767 : END DO
3768 : END IF
3769 167898 : IF (.NOT. found) THEN
3770 21032 : nsize = nsize + 1
3771 21032 : CALL reallocate(array, 1, nsize)
3772 12 : SELECT CASE (fmt)
3773 : CASE (1)
3774 12 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1
3775 : CASE (2)
3776 1489 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2
3777 : CASE (3)
3778 11676 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2, atm3
3779 : CASE (4)
3780 21032 : WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2, atm3, atm4
3781 : END SELECT
3782 : END IF
3783 :
3784 171120 : END SUBROUTINE store_FF_missing_par
3785 :
3786 : ! **************************************************************************************************
3787 : !> \brief Search sorted 2d array of integers for a first occurence of value `val` in row `row`
3788 : !> \param array 2d array of integers
3789 : !> \param val value to search
3790 : !> \param row row to search, default = 1
3791 : !> \return column index if `val` is found in the row `row` of `array`; zero otherwise
3792 : ! **************************************************************************************************
3793 45098 : FUNCTION bsearch_leftmost_2d(array, val, row) RESULT(res)
3794 : INTEGER, INTENT(IN) :: array(:, :), val
3795 : INTEGER, INTENT(IN), OPTIONAL :: row
3796 : INTEGER :: res
3797 :
3798 : INTEGER :: left, locRow, mid, right
3799 :
3800 45098 : locRow = 1
3801 45098 : IF (PRESENT(row)) locRow = row
3802 :
3803 45098 : left = 1
3804 45098 : right = UBOUND(array, dim=2)
3805 :
3806 571050 : DO WHILE (left < right)
3807 525952 : mid = (left + right)/2
3808 571050 : IF (array(locRow, mid) < val) THEN
3809 349610 : left = mid + 1
3810 : ELSE
3811 : right = mid
3812 : END IF
3813 : END DO
3814 :
3815 45098 : res = left
3816 :
3817 : ! Not found:
3818 45098 : IF (array(locRow, res) /= val) res = 0
3819 :
3820 45098 : END FUNCTION bsearch_leftmost_2d
3821 :
3822 : END MODULE force_fields_all
|
