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