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