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