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