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