(git:ed6f26b)
Loading...
Searching...
No Matches
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-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \par History
10!> Splitting and cleaning the original force_field_pack - May 2007
11!> Teodoro Laino - Zurich University
12!> \author CJM
13! **************************************************************************************************
14MODULE 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
116CONTAINS
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) THEN
3297 cpwarn("Multiple "//trim(tag_label)//" declarations: "//trim(item)//" overwriting!")
3298 END IF
3299
3300 END SUBROUTINE issue_duplications
3301
3302! **************************************************************************************************
3303!> \brief Store informations on possible missing ForceFields parameters
3304!> \param atm1 ...
3305!> \param atm2 ...
3306!> \param atm3 ...
3307!> \param atm4 ...
3308!> \param type_name ...
3309!> \param fatal ...
3310!> \param array ...
3311! **************************************************************************************************
3312 SUBROUTINE store_ff_missing_par(atm1, atm2, atm3, atm4, type_name, fatal, array)
3313 CHARACTER(LEN=*), INTENT(IN) :: atm1
3314 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: atm2, atm3, atm4
3315 CHARACTER(LEN=*), INTENT(IN) :: type_name
3316 LOGICAL, INTENT(INOUT), OPTIONAL :: fatal
3317 CHARACTER(LEN=default_string_length), &
3318 DIMENSION(:), POINTER :: array
3319
3320 CHARACTER(LEN=10) :: sfmt
3321 CHARACTER(LEN=9) :: my_atm1, my_atm2, my_atm3, my_atm4
3322 CHARACTER(LEN=default_path_length) :: my_format
3323 INTEGER :: fmt, i, nsize
3324 LOGICAL :: found
3325
3326 nsize = 0
3327 fmt = 1
3328 my_format = '(T2,"FORCEFIELD| Missing ","'//trim(type_name)// &
3329 '",T40,"(",A9,")")'
3330 IF (PRESENT(atm2)) fmt = fmt + 1
3331 IF (PRESENT(atm3)) fmt = fmt + 1
3332 IF (PRESENT(atm4)) fmt = fmt + 1
3333 CALL integer_to_string(fmt - 1, sfmt)
3334 IF (fmt > 1) &
3335 my_format = '(T2,"FORCEFIELD| Missing ","'//trim(type_name)// &
3336 '",T40,"(",A9,'//trim(sfmt)//'(",",A9),")")'
3337 IF (PRESENT(fatal)) fatal = .true.
3338 ! Check for previous already stored equal force fields
3339 IF (ASSOCIATED(array)) nsize = SIZE(array)
3340 found = .false.
3341 IF (nsize >= 1) THEN
3342 DO i = 1, nsize
3343 SELECT CASE (type_name)
3344 CASE ("Bond")
3345 IF (index(array(i) (21:39), "Bond") == 0) cycle
3346 my_atm1 = array(i) (41:49)
3347 my_atm2 = array(i) (51:59)
3348 CALL compress(my_atm1, .true.)
3349 CALL compress(my_atm2, .true.)
3350 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. &
3351 ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .true.
3352 CASE ("Angle")
3353 IF (index(array(i) (21:39), "Angle") == 0) cycle
3354 my_atm1 = array(i) (41:49)
3355 my_atm2 = array(i) (51:59)
3356 my_atm3 = array(i) (61:69)
3357 CALL compress(my_atm1, .true.)
3358 CALL compress(my_atm2, .true.)
3359 CALL compress(my_atm3, .true.)
3360 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. &
3361 ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) &
3362 found = .true.
3363 CASE ("Urey-Bradley")
3364 IF (index(array(i) (21:39), "Urey-Bradley") == 0) cycle
3365 my_atm1 = array(i) (41:49)
3366 my_atm2 = array(i) (51:59)
3367 my_atm3 = array(i) (61:69)
3368 CALL compress(my_atm1, .true.)
3369 CALL compress(my_atm2, .true.)
3370 CALL compress(my_atm3, .true.)
3371 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. &
3372 ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) &
3373 found = .true.
3374 CASE ("Torsion")
3375 IF (index(array(i) (21:39), "Torsion") == 0) cycle
3376 my_atm1 = array(i) (41:49)
3377 my_atm2 = array(i) (51:59)
3378 my_atm3 = array(i) (61:69)
3379 my_atm4 = array(i) (71:79)
3380 CALL compress(my_atm1, .true.)
3381 CALL compress(my_atm2, .true.)
3382 CALL compress(my_atm3, .true.)
3383 CALL compress(my_atm4, .true.)
3384 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3385 ((atm1 == my_atm4) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm1))) &
3386 found = .true.
3387 CASE ("Improper")
3388 IF (index(array(i) (21:39), "Improper") == 0) cycle
3389 my_atm1 = array(i) (41:49)
3390 my_atm2 = array(i) (51:59)
3391 my_atm3 = array(i) (61:69)
3392 my_atm4 = array(i) (71:79)
3393 CALL compress(my_atm1, .true.)
3394 CALL compress(my_atm2, .true.)
3395 CALL compress(my_atm3, .true.)
3396 CALL compress(my_atm4, .true.)
3397 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3398 ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4)) .OR. &
3399 ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3)) .OR. &
3400 ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm2)) .OR. &
3401 ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm3)) .OR. &
3402 ((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3))) &
3403 found = .true.
3404
3405 CASE ("Out of plane bend")
3406 IF (index(array(i) (21:39), "Out of plane bend") == 0) cycle
3407 my_atm1 = array(i) (41:49)
3408 my_atm2 = array(i) (51:59)
3409 my_atm3 = array(i) (61:69)
3410 my_atm4 = array(i) (71:79)
3411 CALL compress(my_atm1, .true.)
3412 CALL compress(my_atm2, .true.)
3413 CALL compress(my_atm3, .true.)
3414 CALL compress(my_atm4, .true.)
3415 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. &
3416 ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4))) &
3417 found = .true.
3418
3419 CASE ("Charge")
3420 IF (index(array(i) (21:39), "Charge") == 0) cycle
3421 my_atm1 = array(i) (41:49)
3422 CALL compress(my_atm1, .true.)
3423 IF (atm1 == my_atm1) found = .true.
3424 CASE ("Spline_Bond_Env", "Spline_Non_Bond_Env")
3425 IF (index(array(i) (21:39), "Spline_") == 0) cycle
3426 fmt = 0
3427 my_atm1 = array(i) (41:49)
3428 my_atm2 = array(i) (51:59)
3429 CALL compress(my_atm1, .true.)
3430 CALL compress(my_atm2, .true.)
3431 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. &
3432 ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .true.
3433 CASE DEFAULT
3434 ! Should never reach this point
3435 cpabort("")
3436 END SELECT
3437 IF (found) EXIT
3438 END DO
3439 END IF
3440 IF (.NOT. found) THEN
3441 nsize = nsize + 1
3442 CALL reallocate(array, 1, nsize)
3443 SELECT CASE (fmt)
3444 CASE (1)
3445 WRITE (array(nsize), fmt=trim(my_format)) atm1
3446 CASE (2)
3447 WRITE (array(nsize), fmt=trim(my_format)) atm1, atm2
3448 CASE (3)
3449 WRITE (array(nsize), fmt=trim(my_format)) atm1, atm2, atm3
3450 CASE (4)
3451 WRITE (array(nsize), fmt=trim(my_format)) atm1, atm2, atm3, atm4
3452 END SELECT
3453 END IF
3454
3455 END SUBROUTINE store_ff_missing_par
3456
3457! **************************************************************************************************
3458!> \brief Search sorted 2d array of integers for a first occurence of value `val` in row `row`
3459!> \param array 2d array of integers
3460!> \param val value to search
3461!> \param row row to search, default = 1
3462!> \return column index if `val` is found in the row `row` of `array`; zero otherwise
3463! **************************************************************************************************
3464 FUNCTION bsearch_leftmost_2d(array, val, row) RESULT(res)
3465 INTEGER, INTENT(IN) :: array(:, :), val
3466 INTEGER, INTENT(IN), OPTIONAL :: row
3467 INTEGER :: res
3468
3469 INTEGER :: left, locrow, mid, right
3470
3471 locrow = 1
3472 IF (PRESENT(row)) locrow = row
3473
3474 left = 1
3475 right = ubound(array, dim=2)
3476
3477 DO WHILE (left < right)
3478 mid = (left + right)/2
3479 IF (array(locrow, mid) < val) THEN
3480 left = mid + 1
3481 ELSE
3482 right = mid
3483 END IF
3484 END DO
3485
3486 res = left
3487
3488 ! Not found:
3489 IF (array(locrow, res) /= val) res = 0
3490
3491 END FUNCTION bsearch_leftmost_2d
3492
3493END MODULE force_fields_all
3494
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
external_potential_types::get_potential
Definition external_potential_types.F:276
external_potential_types::set_potential
Definition external_potential_types.F:296
memory_utilities::reallocate
Definition memory_utilities.F:27
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:181
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:124
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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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
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
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:385
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:1156
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:731
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:704
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:1047
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:704
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:763
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:655
pair_potential_types::pair_potential_lj_create
subroutine, public pair_potential_lj_create(lj)
Cleans the LJ potential type.
Definition pair_potential_types.F:1453
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:942
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:3239
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:3191
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_linked_list_input::cp_sll_val_type
represent a single linked list that stores pointers to the elements
Definition cp_linked_list_input.F:341
damping_dipole_types::damping_p_type
Definition damping_dipole_types.F:45
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:56
external_potential_types::fist_potential_type
Definition external_potential_types.F:85
force_field_kind_types::bend_kind_type
Definition force_field_kind_types.F:58
force_field_kind_types::bond_kind_type
Definition force_field_kind_types.F:51
force_field_kind_types::impr_kind_type
Definition force_field_kind_types.F:83
force_field_kind_types::opbend_kind_type
Definition force_field_kind_types.F:90
force_field_kind_types::torsion_kind_type
Definition force_field_kind_types.F:74
force_field_kind_types::ub_kind_type
Definition force_field_kind_types.F:67
force_field_types::amber_info_type
Definition force_field_types.F:136
force_field_types::charmm_info_type
Definition force_field_types.F:99
force_field_types::force_field_type
Definition force_field_types.F:185
force_field_types::gromos_info_type
Definition force_field_types.F:163
force_field_types::input_info_type
Definition force_field_types.F:38
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
input_val_types::val_type
a type to have a wrapper that stores any basic fortran type
Definition input_val_types.F:63
molecule_kind_types::bend_type
Definition molecule_kind_types.F:76
molecule_kind_types::bond_type
Definition molecule_kind_types.F:69
molecule_kind_types::impr_type
Definition molecule_kind_types.F:97
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_kind_types::opbend_type
Definition molecule_kind_types.F:104
molecule_kind_types::shell_type
Definition molecule_kind_types.F:62
molecule_kind_types::torsion_type
Definition molecule_kind_types.F:90
molecule_kind_types::ub_type
Definition molecule_kind_types.F:83
molecule_types::molecule_type
Definition molecule_types.F:110
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:416
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:383
particle_types::particle_type
Definition particle_types.F:35
qmmm_types_low::qmmm_env_mm_type
...
Definition qmmm_types_low.F:238
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
splines_types::spline_data_p_type
Definition splines_types.F:55
splines_types::spline_environment_type
Definition splines_types.F:65