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