(git:b279b6b)
topology_amber.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Handles all functions used to read and interpret AMBER coordinates
10 !> and topology files
11 !>
12 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
13 ! **************************************************************************************************
14 MODULE topology_amber
15  USE cp_log_handling, ONLY: cp_get_default_logger,&
16  cp_logger_type,&
17  cp_to_string
18  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19  cp_print_key_unit_nr
20  USE cp_parser_methods, ONLY: parser_get_next_line,&
21  parser_get_object,&
22  parser_search_string,&
23  parser_test_next_token
24  USE cp_parser_types, ONLY: cp_parser_type,&
25  parser_create,&
26  parser_release
27  USE cp_units, ONLY: cp_unit_to_cp2k
28  USE force_field_types, ONLY: amber_info_type
29  USE input_cp2k_restarts_util, ONLY: section_velocity_val_set
30  USE input_section_types, ONLY: section_vals_get_subs_vals,&
31  section_vals_type
32  USE kinds, ONLY: default_string_length,&
33  dp
34  USE memory_utilities, ONLY: reallocate
35  USE message_passing, ONLY: mp_para_env_type
36  USE particle_types, ONLY: particle_type
37  USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
38  USE string_table, ONLY: id2str,&
39  s2s,&
40  str2id
41  USE topology_generate_util, ONLY: topology_generate_molname
42  USE topology_types, ONLY: atom_info_type,&
43  connectivity_info_type,&
44  topology_parameters_type
45  USE util, ONLY: sort
46 #include "./base/base_uses.f90"
47 
48  IMPLICIT NONE
49 
50  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_amber'
51  REAL(KIND=dp), PARAMETER, PRIVATE :: amber_conv_factor = 20.4550_dp, &
52  amber_conv_charge = 18.2223_dp
53  INTEGER, PARAMETER, PRIVATE :: buffer_size = 1
54 
55  PRIVATE
56  PUBLIC :: read_coordinate_crd, read_connectivity_amber, rdparm_amber_8
57 
58  ! Reading Amber sections routines
59  INTERFACE rd_amber_section
60  MODULE PROCEDURE rd_amber_section_i1, rd_amber_section_c1, rd_amber_section_r1, &
61  rd_amber_section_i3, rd_amber_section_i4, rd_amber_section_i5
62  END INTERFACE
63 
64 CONTAINS
65 
66 ! **************************************************************************************************
67 !> \brief Reads the `coord' version generated by the PARM or LEaP programs, as
68 !> well as the `restrt' version, resulting from energy minimization or
69 !> molecular dynamics in SANDER or GIBBS. It may contain velocity and
70 !> periodic box information.
71 !>
72 !> Official Format from the AMBER homepage
73 !> FORMAT(20A4) ITITL
74 !> ITITL : the title of the current run, from the AMBER
75 !> parameter/topology file
76 !>
77 !> FORMAT(I5,5E15.7) NATOM,TIME
78 !> NATOM : total number of atoms in coordinate file
79 !> TIME : option, current time in the simulation (picoseconds)
80 !>
81 !> FORMAT(6F12.7) (X(i), Y(i), Z(i), i = 1,NATOM)
82 !> X,Y,Z : coordinates
83 !>
84 !> IF dynamics
85 !>
86 !> FORMAT(6F12.7) (VX(i), VY(i), VZ(i), i = 1,NATOM)
87 !> VX,VY,VZ : velocities (units: Angstroms per 1/20.455 ps)
88 !>
89 !> IF constant pressure (in 4.1, also constant volume)
90 !>
91 !> FORMAT(6F12.7) BOX(1), BOX(2), BOX(3)
92 !> BOX : size of the periodic box
93 !>
94 !>
95 !> \param topology ...
96 !> \param para_env ...
97 !> \param subsys_section ...
98 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
99 ! **************************************************************************************************
100  SUBROUTINE read_coordinate_crd(topology, para_env, subsys_section)
101  TYPE(topology_parameters_type) :: topology
102  TYPE(mp_para_env_type), POINTER :: para_env
103  TYPE(section_vals_type), POINTER :: subsys_section
104 
105  CHARACTER(len=*), PARAMETER :: routinen = 'read_coordinate_crd'
106 
107  CHARACTER(LEN=default_string_length) :: string
108  INTEGER :: handle, iw, j, natom
109  LOGICAL :: my_end, setup_velocities
110  REAL(kind=dp), DIMENSION(:, :), POINTER :: velocity
111  TYPE(atom_info_type), POINTER :: atom_info
112  TYPE(cp_logger_type), POINTER :: logger
113  TYPE(cp_parser_type) :: parser
114  TYPE(section_vals_type), POINTER :: velocity_section
115 
116  NULLIFY (logger, velocity)
117  logger => cp_get_default_logger()
118  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/CRD_INFO", &
119  extension=".subsysLog")
120  CALL timeset(routinen, handle)
121 
122  atom_info => topology%atom_info
123  IF (iw > 0) WRITE (iw, *) " Reading in CRD file ", trim(topology%coord_file_name)
124 
125  ! Title Section
126  IF (iw > 0) WRITE (iw, '(T2,A)') 'CRD_INFO| Parsing the TITLE section'
127  CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
128  CALL parser_get_next_line(parser, 1)
129  ! Title may be missing
130  IF (parser_test_next_token(parser) == "STR") THEN
131  CALL parser_get_object(parser, string, string_length=default_string_length)
132  IF (iw > 0) WRITE (iw, '(T2,A)') 'CRD_INFO| '//trim(string)
133  ! Natom and Time (which we ignore)
134  CALL parser_get_next_line(parser, 1)
135  END IF
136  CALL parser_get_object(parser, natom)
137  topology%natoms = natom
138  IF (iw > 0) WRITE (iw, '(T2,A,I0)') 'CRD_INFO| Number of atoms: ', natom
139  CALL reallocate(atom_info%id_molname, 1, natom)
140  CALL reallocate(atom_info%id_resname, 1, natom)
141  CALL reallocate(atom_info%resid, 1, natom)
142  CALL reallocate(atom_info%id_atmname, 1, natom)
143  CALL reallocate(atom_info%r, 1, 3, 1, natom)
144  CALL reallocate(atom_info%atm_mass, 1, natom)
145  CALL reallocate(atom_info%atm_charge, 1, natom)
146  CALL reallocate(atom_info%occup, 1, natom)
147  CALL reallocate(atom_info%beta, 1, natom)
148  CALL reallocate(atom_info%id_element, 1, natom)
149 
150  ! Element is assigned on the basis of the atm_name
151  topology%aa_element = .true.
152 
153  ! Coordinates
154  CALL parser_get_next_line(parser, 1, at_end=my_end)
155  DO j = 1, natom - mod(natom, 2), 2
156  IF (my_end) EXIT
157  READ (parser%input_line, *) atom_info%r(1, j), atom_info%r(2, j), atom_info%r(3, j), &
158  atom_info%r(1, j + 1), atom_info%r(2, j + 1), atom_info%r(3, j + 1)
159  ! All these information will have to be setup elsewhere..
160  ! CRD file does not contain anything related..
161  atom_info%id_atmname(j) = str2id(s2s("__UNDEF__"))
162  atom_info%id_molname(j) = str2id(s2s("__UNDEF__"))
163  atom_info%id_resname(j) = str2id(s2s("__UNDEF__"))
164  atom_info%id_element(j) = str2id(s2s("__UNDEF__"))
165  atom_info%resid(j) = huge(0)
166  atom_info%atm_mass(j) = huge(0.0_dp)
167  atom_info%atm_charge(j) = -huge(0.0_dp)
168  atom_info%r(1, j) = cp_unit_to_cp2k(atom_info%r(1, j), "angstrom")
169  atom_info%r(2, j) = cp_unit_to_cp2k(atom_info%r(2, j), "angstrom")
170  atom_info%r(3, j) = cp_unit_to_cp2k(atom_info%r(3, j), "angstrom")
171 
172  atom_info%id_atmname(j + 1) = str2id(s2s("__UNDEF__"))
173  atom_info%id_molname(j + 1) = str2id(s2s("__UNDEF__"))
174  atom_info%id_resname(j + 1) = str2id(s2s("__UNDEF__"))
175  atom_info%id_element(j + 1) = str2id(s2s("__UNDEF__"))
176  atom_info%resid(j + 1) = huge(0)
177  atom_info%atm_mass(j + 1) = huge(0.0_dp)
178  atom_info%atm_charge(j + 1) = -huge(0.0_dp)
179  atom_info%r(1, j + 1) = cp_unit_to_cp2k(atom_info%r(1, j + 1), "angstrom")
180  atom_info%r(2, j + 1) = cp_unit_to_cp2k(atom_info%r(2, j + 1), "angstrom")
181  atom_info%r(3, j + 1) = cp_unit_to_cp2k(atom_info%r(3, j + 1), "angstrom")
182 
183  CALL parser_get_next_line(parser, 1, at_end=my_end)
184  END DO
185  ! Trigger error
186  IF ((my_end) .AND. (j /= natom - mod(natom, 2) + 1)) THEN
187  IF (j /= natom) &
188  cpabort("Error while reading CRD file. Unexpected end of file.")
189  ELSE IF (mod(natom, 2) /= 0) THEN
190  ! In case let's handle the last atom
191  j = natom
192  READ (parser%input_line, *) atom_info%r(1, j), atom_info%r(2, j), atom_info%r(3, j)
193  ! All these information will have to be setup elsewhere..
194  ! CRD file does not contain anything related..
195  atom_info%id_atmname(j) = str2id(s2s("__UNDEF__"))
196  atom_info%id_molname(j) = str2id(s2s("__UNDEF__"))
197  atom_info%id_resname(j) = str2id(s2s("__UNDEF__"))
198  atom_info%id_element(j) = str2id(s2s("__UNDEF__"))
199  atom_info%resid(j) = huge(0)
200  atom_info%atm_mass(j) = huge(0.0_dp)
201  atom_info%atm_charge(j) = -huge(0.0_dp)
202  atom_info%r(1, j) = cp_unit_to_cp2k(atom_info%r(1, j), "angstrom")
203  atom_info%r(2, j) = cp_unit_to_cp2k(atom_info%r(2, j), "angstrom")
204  atom_info%r(3, j) = cp_unit_to_cp2k(atom_info%r(3, j), "angstrom")
205 
206  CALL parser_get_next_line(parser, 1, at_end=my_end)
207  END IF
208 
209  IF (my_end) THEN
210  IF (j /= natom) &
211  cpwarn("No VELOCITY or BOX information found in CRD file. ")
212  ELSE
213  ! Velocities
214  CALL reallocate(velocity, 1, 3, 1, natom)
215  DO j = 1, natom - mod(natom, 2), 2
216  IF (my_end) EXIT
217  READ (parser%input_line, *) velocity(1, j), velocity(2, j), velocity(3, j), &
218  velocity(1, j + 1), velocity(2, j + 1), velocity(3, j + 1)
219 
220  velocity(1, j) = cp_unit_to_cp2k(velocity(1, j), "angstrom*ps^-1")
221  velocity(2, j) = cp_unit_to_cp2k(velocity(2, j), "angstrom*ps^-1")
222  velocity(3, j) = cp_unit_to_cp2k(velocity(3, j), "angstrom*ps^-1")
223  velocity(1:3, j) = velocity(1:3, j)*amber_conv_factor
224 
225  velocity(1, j + 1) = cp_unit_to_cp2k(velocity(1, j + 1), "angstrom*ps^-1")
226  velocity(2, j + 1) = cp_unit_to_cp2k(velocity(2, j + 1), "angstrom*ps^-1")
227  velocity(3, j + 1) = cp_unit_to_cp2k(velocity(3, j + 1), "angstrom*ps^-1")
228  velocity(1:3, j + 1) = velocity(1:3, j + 1)*amber_conv_factor
229 
230  CALL parser_get_next_line(parser, 1, at_end=my_end)
231  END DO
232  setup_velocities = .true.
233  IF ((my_end) .AND. (j /= natom - mod(natom, 2) + 1)) THEN
234  IF (j /= natom) &
235  CALL cp_warn(__location__, &
236  "No VELOCITY information found in CRD file. Ignoring BOX information. "// &
237  "Please provide the BOX information directly from the main CP2K input! ")
238  setup_velocities = .false.
239  ELSE IF (mod(natom, 2) /= 0) THEN
240  ! In case let's handle the last atom
241  j = natom
242  READ (parser%input_line, *) velocity(1, j), velocity(2, j), velocity(3, j)
243 
244  velocity(1, j) = cp_unit_to_cp2k(velocity(1, j), "angstrom*ps^-1")
245  velocity(2, j) = cp_unit_to_cp2k(velocity(2, j), "angstrom*ps^-1")
246  velocity(3, j) = cp_unit_to_cp2k(velocity(3, j), "angstrom*ps^-1")
247  velocity(1:3, j) = velocity(1:3, j)*amber_conv_factor
248 
249  CALL parser_get_next_line(parser, 1, at_end=my_end)
250  END IF
251  IF (setup_velocities) THEN
252  velocity_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
253  CALL section_velocity_val_set(velocity_section, velocity=velocity, &
254  conv_factor=1.0_dp)
255  END IF
256  DEALLOCATE (velocity)
257  END IF
258  IF (my_end) THEN
259  IF (j /= natom) &
260  cpwarn("BOX information missing in CRD file. ")
261  ELSE
262  IF (j /= natom) &
263  CALL cp_warn(__location__, &
264  "BOX information found in CRD file. They will be ignored. "// &
265  "Please provide the BOX information directly from the main CP2K input!")
266  END IF
267  CALL parser_release(parser)
268  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
269  "PRINT%TOPOLOGY_INFO/CRD_INFO")
270  CALL timestop(handle)
271 
272  END SUBROUTINE read_coordinate_crd
273 
274 ! **************************************************************************************************
275 !> \brief Read AMBER topology file (.top) : At this level we parse only the
276 !> connectivity info the .top file. ForceField information will be
277 !> handled later
278 !>
279 !> \param filename ...
280 !> \param topology ...
281 !> \param para_env ...
282 !> \param subsys_section ...
283 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
284 ! **************************************************************************************************
285  SUBROUTINE read_connectivity_amber(filename, topology, para_env, subsys_section)
286  CHARACTER(LEN=*), INTENT(IN) :: filename
287  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
288  TYPE(mp_para_env_type), POINTER :: para_env
289  TYPE(section_vals_type), POINTER :: subsys_section
290 
291  CHARACTER(len=*), PARAMETER :: routinen = 'read_connectivity_amber'
292 
293  INTEGER :: handle, iw
294  TYPE(atom_info_type), POINTER :: atom_info
295  TYPE(connectivity_info_type), POINTER :: conn_info
296  TYPE(cp_logger_type), POINTER :: logger
297 
298  NULLIFY (logger)
299  CALL timeset(routinen, handle)
300  logger => cp_get_default_logger()
301  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/AMBER_INFO", &
302  extension=".subsysLog")
303 
304  atom_info => topology%atom_info
305  conn_info => topology%conn_info
306 
307  ! Read the Amber topology file
308  CALL rdparm_amber_8(filename, iw, para_env, do_connectivity=.true., do_forcefield=.false., &
309  atom_info=atom_info, conn_info=conn_info)
310 
311  ! Molnames have been internally generated
312  topology%molname_generated = .true.
313 
314  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
315  "PRINT%TOPOLOGY_INFO/AMBER_INFO")
316  CALL timestop(handle)
317  END SUBROUTINE read_connectivity_amber
318 
319 ! **************************************************************************************************
320 !> \brief Access information form the AMBER topology file
321 !> Notes on file structure:
322 !>
323 !> NATOM ! Total number of Atoms
324 !> NTYPES ! Total number of distinct atom types
325 !> NBONH ! Number of bonds containing hydrogens
326 !> MBONA ! Number of bonds not containing hydrogens
327 !> NTHETH ! Number of angles containing hydrogens
328 !> MTHETA ! Number of angles not containing hydrogens
329 !> NPHIH ! Number of dihedrals containing hydrogens
330 !> MPHIA ! Number of dihedrals not containing hydrogens
331 !> NHPARM ! currently NOT USED
332 !> NPARM ! set to 1 if LES is used
333 !> NNB ! number of excluded atoms
334 !> NRES ! Number of residues
335 !> NBONA ! MBONA + number of constraint bonds ( in v.8 NBONA=MBONA)
336 !> NTHETA ! MTHETA + number of constraint angles ( in v.8 NBONA=MBONA)
337 !> NPHIA ! MPHIA + number of constraint dihedrals ( in v.8 NBONA=MBONA)
338 !> NUMBND ! Number of unique bond types
339 !> NUMANG ! Number of unique angle types
340 !> NPTRA ! Number of unique dihedral types
341 !> NATYP ! Number of atom types in parameter file
342 !> NPHB ! Number of distinct 10-12 hydrogen bond pair types
343 !> IFPERT ! Variable not used in this converter...
344 !> NBPER ! Variable not used in this converter...
345 !> NGPER ! Variable not used in this converter...
346 !> NDPER ! Variable not used in this converter...
347 !> MBPER ! Variable not used in this converter...
348 !> MGPER ! Variable not used in this converter...
349 !> MDPER ! Variable not used in this converter...
350 !> IFBOX ! Variable not used in this converter...
351 !> NMXRS ! Variable not used in this converter...
352 !> IFCAP ! Variable not used in this converter...
353 !> NUMEXTRA ! Variable not used in this converter...
354 !>
355 !> \param filename ...
356 !> \param output_unit ...
357 !> \param para_env ...
358 !> \param do_connectivity ...
359 !> \param do_forcefield ...
360 !> \param atom_info ...
361 !> \param conn_info ...
362 !> \param amb_info ...
363 !> \param particle_set ...
364 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
365 ! **************************************************************************************************
366  SUBROUTINE rdparm_amber_8(filename, output_unit, para_env, do_connectivity, &
367  do_forcefield, atom_info, conn_info, amb_info, particle_set)
368 
369  CHARACTER(LEN=*), INTENT(IN) :: filename
370  INTEGER, INTENT(IN) :: output_unit
371  TYPE(mp_para_env_type), POINTER :: para_env
372  LOGICAL, INTENT(IN) :: do_connectivity, do_forcefield
373  TYPE(atom_info_type), OPTIONAL, POINTER :: atom_info
374  TYPE(connectivity_info_type), OPTIONAL, POINTER :: conn_info
375  TYPE(amber_info_type), OPTIONAL, POINTER :: amb_info
376  TYPE(particle_type), DIMENSION(:), OPTIONAL, &
377  POINTER :: particle_set
378 
379  CHARACTER(len=*), PARAMETER :: routinen = 'rdparm_amber_8'
380 
381  CHARACTER(LEN=default_string_length) :: input_format, section
382  CHARACTER(LEN=default_string_length), &
383  ALLOCATABLE, DIMENSION(:) :: isymbl, labres, strtmp_a
384  INTEGER :: handle, handle2, i, ifbox, ifcap, ifpert, index_now, info(31), istart, mbona, &
385  mbper, mdper, mgper, mphia, mtheta, natom, natom_prev, natyp, nbona, nbond_prev, nbonh, &
386  nbper, ndper, ngper, nhparm, nmxrs, nnb, nparm, nphb, nphi_prev, nphia, nphih, nptra, &
387  nres, nsize, ntheta, ntheta_prev, ntheth, ntypes, numang, numbnd, numextra, &
388  unique_torsions
389  INTEGER, ALLOCATABLE, DIMENSION(:) :: iac, ib, ibh, icb, icbh, ico, icp, icph, &
390  ict, icth, ip, iph, ipres, it, ith, &
391  iwork, jb, jbh, jp, jph, jt, jth, kp, &
392  kph, kt, kth, lp, lph
393  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: full_torsions
394  LOGICAL :: check, valid_format
395  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: asol, bsol, cn1, cn2, phase, pk, pn, &
396  req, rk, teq, tk
397  TYPE(cp_parser_type) :: parser
398 
399  CALL timeset(routinen, handle)
400  IF (output_unit > 0) WRITE (output_unit, '(/,A)') " AMBER_INFO| Reading Amber Topology File: "// &
401  trim(filename)
402  CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.true.)
403  valid_format = check_amber_8_std(parser, output_unit)
404  IF (valid_format) THEN
405  DO WHILE (get_section_parmtop(parser, section, input_format))
406  SELECT CASE (trim(section))
407  CASE ("TITLE")
408  ! Who cares about the title?
409  cycle
410  CASE ("POINTERS")
411  CALL rd_amber_section(parser, section, info, 31)
412  ! Assign pointers to the corresponding labels
413  ! just for convenience to have something more human readable
414  natom = info(1)
415  ntypes = info(2)
416  nbonh = info(3)
417  mbona = info(4)
418  ntheth = info(5)
419  mtheta = info(6)
420  nphih = info(7)
421  mphia = info(8)
422  nhparm = info(9)
423  nparm = info(10)
424  nnb = info(11)
425  nres = info(12)
426  nbona = info(13)
427  ntheta = info(14)
428  nphia = info(15)
429  numbnd = info(16)
430  numang = info(17)
431  nptra = info(18)
432  natyp = info(19)
433  nphb = info(20)
434  ifpert = info(21)
435  nbper = info(22)
436  ngper = info(23)
437  ndper = info(24)
438  mbper = info(25)
439  mgper = info(26)
440  mdper = info(27)
441  ifbox = info(28)
442  nmxrs = info(29)
443  ifcap = info(30)
444  numextra = info(31)
445 
446  ! Print some info if requested
447  IF (output_unit > 0) THEN
448  WRITE (output_unit, '(A,/)') " AMBER_INFO| Information from AMBER topology file:"
449  WRITE (output_unit, 1000) &
450  natom, ntypes, nbonh, mbona, ntheth, mtheta, nphih, &
451  mphia, nhparm, nparm, nnb, nres, nbona, ntheta, &
452  nphia, numbnd, numang, nptra, natyp, nphb, ifbox, &
453  nmxrs, ifcap, numextra
454  END IF
455 
456  ! Allocate temporary arrays
457  IF (do_connectivity) THEN
458  check = PRESENT(atom_info) .AND. PRESENT(conn_info)
459  cpassert(check)
460  natom_prev = 0
461  IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname)
462  ! Allocate for extracting connectivity infos
463  ALLOCATE (labres(nres))
464  ALLOCATE (ipres(nres))
465  END IF
466  IF (do_forcefield) THEN
467  ! Allocate for extracting forcefield infos
468  ALLOCATE (iac(natom))
469  ALLOCATE (ico(ntypes*ntypes))
470  ALLOCATE (rk(numbnd))
471  ALLOCATE (req(numbnd))
472  ALLOCATE (tk(numang))
473  ALLOCATE (teq(numang))
474  ALLOCATE (pk(nptra))
475  ALLOCATE (pn(nptra))
476  ALLOCATE (phase(nptra))
477  ALLOCATE (cn1(ntypes*(ntypes + 1)/2))
478  ALLOCATE (cn2(ntypes*(ntypes + 1)/2))
479  ALLOCATE (asol(ntypes*(ntypes + 1)/2))
480  ALLOCATE (bsol(ntypes*(ntypes + 1)/2))
481  END IF
482  ! Always Allocate
483  ALLOCATE (ibh(nbonh))
484  ALLOCATE (jbh(nbonh))
485  ALLOCATE (icbh(nbonh))
486  ALLOCATE (ib(nbona))
487  ALLOCATE (jb(nbona))
488  ALLOCATE (icb(nbona))
489  ALLOCATE (ith(ntheth))
490  ALLOCATE (jth(ntheth))
491  ALLOCATE (kth(ntheth))
492  ALLOCATE (icth(ntheth))
493  ALLOCATE (it(ntheta))
494  ALLOCATE (jt(ntheta))
495  ALLOCATE (kt(ntheta))
496  ALLOCATE (ict(ntheta))
497  ALLOCATE (iph(nphih))
498  ALLOCATE (jph(nphih))
499  ALLOCATE (kph(nphih))
500  ALLOCATE (lph(nphih))
501  ALLOCATE (icph(nphih))
502  ALLOCATE (ip(nphia))
503  ALLOCATE (jp(nphia))
504  ALLOCATE (kp(nphia))
505  ALLOCATE (lp(nphia))
506  ALLOCATE (icp(nphia))
507  CASE ("ATOM_NAME")
508  ! Atom names are just ignored according the CP2K philosophy
509  cycle
510  CASE ("AMBER_ATOM_TYPE")
511  IF (.NOT. do_connectivity) cycle
512  CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom)
513  ALLOCATE (strtmp_a(natom))
514  CALL rd_amber_section(parser, section, strtmp_a, natom)
515  DO i = 1, natom
516  atom_info%id_atmname(natom_prev + i) = str2id(strtmp_a(i))
517  END DO
518  DEALLOCATE (strtmp_a)
519  CASE ("CHARGE")
520  IF (.NOT. do_connectivity) cycle
521  CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom)
522  CALL rd_amber_section(parser, section, atom_info%atm_charge(natom_prev + 1:), natom)
523  ! Convert charges into atomic units
524  atom_info%atm_charge(natom_prev + 1:) = atom_info%atm_charge(natom_prev + 1:)/amber_conv_charge
525  CASE ("MASS")
526  IF (.NOT. do_connectivity) cycle
527  CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom)
528  CALL rd_amber_section(parser, section, atom_info%atm_mass(natom_prev + 1:), natom)
529  CASE ("RESIDUE_LABEL")
530  IF (.NOT. do_connectivity) cycle
531  CALL reallocate(atom_info%id_resname, 1, natom_prev + natom)
532  CALL rd_amber_section(parser, section, labres, nres)
533  CASE ("RESIDUE_POINTER")
534  IF (.NOT. do_connectivity) cycle
535  CALL reallocate(atom_info%resid, 1, natom_prev + natom)
536  CALL rd_amber_section(parser, section, ipres, nres)
537  CASE ("ATOM_TYPE_INDEX")
538  IF (.NOT. do_forcefield) cycle
539  CALL rd_amber_section(parser, section, iac, natom)
540  CASE ("NONBONDED_PARM_INDEX")
541  IF (.NOT. do_forcefield) cycle
542  CALL rd_amber_section(parser, section, ico, ntypes**2)
543  CASE ("BOND_FORCE_CONSTANT")
544  IF (.NOT. do_forcefield) cycle
545  CALL rd_amber_section(parser, section, rk, numbnd)
546  CASE ("BOND_EQUIL_VALUE")
547  IF (.NOT. do_forcefield) cycle
548  CALL rd_amber_section(parser, section, req, numbnd)
549  CASE ("ANGLE_FORCE_CONSTANT")
550  IF (.NOT. do_forcefield) cycle
551  CALL rd_amber_section(parser, section, tk, numang)
552  CASE ("ANGLE_EQUIL_VALUE")
553  IF (.NOT. do_forcefield) cycle
554  CALL rd_amber_section(parser, section, teq, numang)
555  CASE ("DIHEDRAL_FORCE_CONSTANT")
556  IF (.NOT. do_forcefield) cycle
557  CALL rd_amber_section(parser, section, pk, nptra)
558  IF (nptra <= 0) cycle
559  ! Save raw values
560  IF (ASSOCIATED(amb_info%raw_torsion_k)) DEALLOCATE (amb_info%raw_torsion_k)
561  ALLOCATE (amb_info%raw_torsion_k(nptra), source=pk)
562  CASE ("DIHEDRAL_PERIODICITY")
563  IF (.NOT. do_forcefield) cycle
564  CALL rd_amber_section(parser, section, pn, nptra)
565  IF (nptra <= 0) cycle
566  ! Save raw values
567  IF (ASSOCIATED(amb_info%raw_torsion_m)) DEALLOCATE (amb_info%raw_torsion_m)
568  ALLOCATE (amb_info%raw_torsion_m(nptra), source=pn)
569  CASE ("DIHEDRAL_PHASE")
570  IF (.NOT. do_forcefield) cycle
571  CALL rd_amber_section(parser, section, phase, nptra)
572  IF (nptra <= 0) cycle
573  ! Save raw values
574  IF (ASSOCIATED(amb_info%raw_torsion_phi0)) DEALLOCATE (amb_info%raw_torsion_phi0)
575  ALLOCATE (amb_info%raw_torsion_phi0(nptra), source=phase)
576  CASE ("LENNARD_JONES_ACOEF")
577  IF (.NOT. do_forcefield) cycle
578  CALL rd_amber_section(parser, section, cn1, ntypes*(ntypes + 1)/2)
579  CASE ("LENNARD_JONES_BCOEF")
580  IF (.NOT. do_forcefield) cycle
581  CALL rd_amber_section(parser, section, cn2, ntypes*(ntypes + 1)/2)
582  CASE ("HBOND_ACOEF")
583  IF (.NOT. do_forcefield) cycle
584  CALL rd_amber_section(parser, section, asol, nphb)
585  CASE ("HBOND_BCOEF")
586  IF (.NOT. do_forcefield) cycle
587  CALL rd_amber_section(parser, section, bsol, nphb)
588  CASE ("BONDS_INC_HYDROGEN")
589  ! We always need to parse this information both for connectivity and forcefields
590  CALL rd_amber_section(parser, section, ibh, jbh, icbh, nbonh)
591  ! Conver to an atomic index
592  ibh(:) = ibh(:)/3 + 1
593  jbh(:) = jbh(:)/3 + 1
594  CASE ("BONDS_WITHOUT_HYDROGEN")
595  ! We always need to parse this information both for connectivity and forcefields
596  CALL rd_amber_section(parser, section, ib, jb, icb, nbona)
597  ! Conver to an atomic index
598  ib(:) = ib(:)/3 + 1
599  jb(:) = jb(:)/3 + 1
600  CASE ("ANGLES_INC_HYDROGEN")
601  ! We always need to parse this information both for connectivity and forcefields
602  CALL rd_amber_section(parser, section, ith, jth, kth, icth, ntheth)
603  ! Conver to an atomic index
604  ith(:) = ith(:)/3 + 1
605  jth(:) = jth(:)/3 + 1
606  kth(:) = kth(:)/3 + 1
607  CASE ("ANGLES_WITHOUT_HYDROGEN")
608  ! We always need to parse this information both for connectivity and forcefields
609  CALL rd_amber_section(parser, section, it, jt, kt, ict, ntheta)
610  ! Conver to an atomic index
611  it(:) = it(:)/3 + 1
612  jt(:) = jt(:)/3 + 1
613  kt(:) = kt(:)/3 + 1
614  CASE ("DIHEDRALS_INC_HYDROGEN")
615  ! We always need to parse this information both for connectivity and forcefields
616  CALL rd_amber_section(parser, section, iph, jph, kph, lph, icph, nphih)
617  ! Conver to an atomic index
618  iph(:) = iph(:)/3 + 1
619  jph(:) = jph(:)/3 + 1
620  kph(:) = abs(kph(:))/3 + 1
621  lph(:) = abs(lph(:))/3 + 1
622  CASE ("DIHEDRALS_WITHOUT_HYDROGEN")
623  ! We always need to parse this information both for connectivity and forcefields
624  CALL rd_amber_section(parser, section, ip, jp, kp, lp, icp, nphia)
625  ! Conver to an atomic index
626  ip(:) = ip(:)/3 + 1
627  jp(:) = jp(:)/3 + 1
628  kp(:) = abs(kp(:))/3 + 1
629  lp(:) = abs(lp(:))/3 + 1
630  CASE DEFAULT
631  ! Just Ignore other sections...
632  END SELECT
633  END DO
634  ! Save raw torsion info: atom indices and dihedral index
635  IF (do_forcefield .AND. (nphih + nphia > 0)) THEN
636  IF (ASSOCIATED(amb_info%raw_torsion_id)) DEALLOCATE (amb_info%raw_torsion_id)
637  ALLOCATE (amb_info%raw_torsion_id(5, nphih + nphia))
638  DO i = 1, nphih
639  amb_info%raw_torsion_id(1, i) = iph(i)
640  amb_info%raw_torsion_id(2, i) = jph(i)
641  amb_info%raw_torsion_id(3, i) = kph(i)
642  amb_info%raw_torsion_id(4, i) = lph(i)
643  amb_info%raw_torsion_id(5, i) = icph(i)
644  END DO
645  DO i = 1, nphia
646  amb_info%raw_torsion_id(1, nphih + i) = ip(i)
647  amb_info%raw_torsion_id(2, nphih + i) = jp(i)
648  amb_info%raw_torsion_id(3, nphih + i) = kp(i)
649  amb_info%raw_torsion_id(4, nphih + i) = lp(i)
650  amb_info%raw_torsion_id(5, nphih + i) = icp(i)
651  END DO
652  END IF
653  END IF
654 
655  ! Extracts connectivity info from the AMBER topology file
656  IF (do_connectivity) THEN
657  CALL timeset(trim(routinen)//"_connectivity", handle2)
658  ! ----------------------------------------------------------
659  ! Conform Amber Names with CHARMM convention (kind<->charge)
660  ! ----------------------------------------------------------
661  ALLOCATE (isymbl(natom))
662  ALLOCATE (iwork(natom))
663 
664  DO i = 1, SIZE(isymbl)
665  isymbl(i) = id2str(atom_info%id_atmname(natom_prev + i))
666  END DO
667 
668  ! Sort atom names + charges and identify unique types
669  CALL sort(isymbl, natom, iwork)
670 
671  istart = 1
672  DO i = 2, natom
673  IF (trim(isymbl(i)) /= trim(isymbl(istart))) THEN
674  CALL conform_atom_type_low(isymbl, iwork, i, istart, atom_info%atm_charge(natom_prev + 1:))
675  istart = i
676  END IF
677  END DO
678  CALL conform_atom_type_low(isymbl, iwork, i, istart, atom_info%atm_charge(natom_prev + 1:))
679 
680  ! Copy back the modified and conformed atom types
681  DO i = 1, natom
682  atom_info%id_atmname(natom_prev + iwork(i)) = str2id(s2s(isymbl(i)))
683  END DO
684 
685  ! -----------------------------------------------------------
686  ! Fill residue_name and residue_id information before exiting
687  ! -----------------------------------------------------------
688  DO i = 1, nres - 1
689  atom_info%id_resname(natom_prev + ipres(i):natom_prev + ipres(i + 1)) = str2id(s2s(labres(i)))
690  atom_info%resid(natom_prev + ipres(i):natom_prev + ipres(i + 1)) = i
691  END DO
692  atom_info%id_resname(natom_prev + ipres(i):natom_prev + natom) = str2id(s2s(labres(i)))
693  atom_info%resid(natom_prev + ipres(i):natom_prev + natom) = i
694 
695  ! Deallocate when extracting connectivity infos
696  DEALLOCATE (iwork)
697  DEALLOCATE (isymbl)
698  DEALLOCATE (labres)
699  DEALLOCATE (ipres)
700 
701  ! ----------------------------------------------------------
702  ! Copy connectivity
703  ! ----------------------------------------------------------
704  ! BONDS
705  nbond_prev = 0
706  IF (ASSOCIATED(conn_info%bond_a)) nbond_prev = SIZE(conn_info%bond_a)
707 
708  CALL reallocate(conn_info%bond_a, 1, nbond_prev + nbonh + nbona)
709  CALL reallocate(conn_info%bond_b, 1, nbond_prev + nbonh + nbona)
710  DO i = 1, nbonh
711  index_now = nbond_prev + i
712  conn_info%bond_a(index_now) = natom_prev + ibh(i)
713  conn_info%bond_b(index_now) = natom_prev + jbh(i)
714  END DO
715  DO i = 1, nbona
716  index_now = nbond_prev + i + nbonh
717  conn_info%bond_a(index_now) = natom_prev + ib(i)
718  conn_info%bond_b(index_now) = natom_prev + jb(i)
719  END DO
720 
721  ! ANGLES
722  ntheta_prev = 0
723  IF (ASSOCIATED(conn_info%theta_a)) ntheta_prev = SIZE(conn_info%theta_a)
724 
725  CALL reallocate(conn_info%theta_a, 1, ntheta_prev + ntheth + ntheta)
726  CALL reallocate(conn_info%theta_b, 1, ntheta_prev + ntheth + ntheta)
727  CALL reallocate(conn_info%theta_c, 1, ntheta_prev + ntheth + ntheta)
728  DO i = 1, ntheth
729  index_now = ntheta_prev + i
730  conn_info%theta_a(index_now) = natom_prev + ith(i)
731  conn_info%theta_b(index_now) = natom_prev + jth(i)
732  conn_info%theta_c(index_now) = natom_prev + kth(i)
733  END DO
734  DO i = 1, ntheta
735  index_now = ntheta_prev + i + ntheth
736  conn_info%theta_a(index_now) = natom_prev + it(i)
737  conn_info%theta_b(index_now) = natom_prev + jt(i)
738  conn_info%theta_c(index_now) = natom_prev + kt(i)
739  END DO
740 
741  ! TORSIONS
742  ! For torsions we need to find out the unique torsions
743  ! defined in the amber parmtop
744  nphi_prev = 0
745  IF (ASSOCIATED(conn_info%phi_a)) nphi_prev = SIZE(conn_info%phi_a)
746 
747  CALL reallocate(conn_info%phi_a, 1, nphi_prev + nphih + nphia)
748  CALL reallocate(conn_info%phi_b, 1, nphi_prev + nphih + nphia)
749  CALL reallocate(conn_info%phi_c, 1, nphi_prev + nphih + nphia)
750  CALL reallocate(conn_info%phi_d, 1, nphi_prev + nphih + nphia)
751 
752  IF (nphih + nphia /= 0) THEN
753  ALLOCATE (full_torsions(4, nphih + nphia))
754  ALLOCATE (iwork(nphih + nphia))
755 
756  DO i = 1, nphih
757  full_torsions(1, i) = iph(i)
758  full_torsions(2, i) = jph(i)
759  full_torsions(3, i) = kph(i)
760  full_torsions(4, i) = lph(i)
761  END DO
762  DO i = 1, nphia
763  full_torsions(1, nphih + i) = ip(i)
764  full_torsions(2, nphih + i) = jp(i)
765  full_torsions(3, nphih + i) = kp(i)
766  full_torsions(4, nphih + i) = lp(i)
767  END DO
768  CALL sort(full_torsions, 1, nphih + nphia, 1, 4, iwork)
769 
770  unique_torsions = nphi_prev + 1
771  conn_info%phi_a(unique_torsions) = natom_prev + full_torsions(1, 1)
772  conn_info%phi_b(unique_torsions) = natom_prev + full_torsions(2, 1)
773  conn_info%phi_c(unique_torsions) = natom_prev + full_torsions(3, 1)
774  conn_info%phi_d(unique_torsions) = natom_prev + full_torsions(4, 1)
775  DO i = 2, nphih + nphia
776  IF ((full_torsions(1, i) /= full_torsions(1, i - 1)) .OR. &
777  (full_torsions(2, i) /= full_torsions(2, i - 1)) .OR. &
778  (full_torsions(3, i) /= full_torsions(3, i - 1)) .OR. &
779  (full_torsions(4, i) /= full_torsions(4, i - 1))) THEN
780  unique_torsions = unique_torsions + 1
781  conn_info%phi_a(unique_torsions) = natom_prev + full_torsions(1, i)
782  conn_info%phi_b(unique_torsions) = natom_prev + full_torsions(2, i)
783  conn_info%phi_c(unique_torsions) = natom_prev + full_torsions(3, i)
784  conn_info%phi_d(unique_torsions) = natom_prev + full_torsions(4, i)
785  END IF
786  END DO
787  CALL reallocate(conn_info%phi_a, 1, unique_torsions)
788  CALL reallocate(conn_info%phi_b, 1, unique_torsions)
789  CALL reallocate(conn_info%phi_c, 1, unique_torsions)
790  CALL reallocate(conn_info%phi_d, 1, unique_torsions)
791 
792  DEALLOCATE (full_torsions)
793  DEALLOCATE (iwork)
794  END IF
795  ! IMPROPERS
796  CALL reallocate(conn_info%impr_a, 1, 0)
797  CALL reallocate(conn_info%impr_b, 1, 0)
798  CALL reallocate(conn_info%impr_c, 1, 0)
799  CALL reallocate(conn_info%impr_d, 1, 0)
800 
801  ! ----------------------------------------------------------
802  ! Generate molecule names
803  ! ----------------------------------------------------------
804  CALL reallocate(atom_info%id_molname, 1, natom_prev + natom)
805  atom_info%id_molname(natom_prev + 1:natom_prev + natom) = str2id(s2s("__UNDEF__"))
806  CALL topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, &
807  atom_info%id_molname(natom_prev + 1:natom_prev + natom))
808  CALL timestop(handle2)
809  END IF
810 
811  ! Extracts force fields info from the AMBER topology file
812  IF (do_forcefield) THEN
813  CALL timeset(trim(routinen)//"_forcefield", handle2)
814  ! ----------------------------------------------------------
815  ! Force Fields informations related to bonds
816  ! ----------------------------------------------------------
817  CALL reallocate(amb_info%bond_a, 1, buffer_size)
818  CALL reallocate(amb_info%bond_b, 1, buffer_size)
819  CALL reallocate(amb_info%bond_k, 1, buffer_size)
820  CALL reallocate(amb_info%bond_r0, 1, buffer_size)
821  nsize = 0
822  ! Bonds containing hydrogens
823  CALL post_process_bonds_info(amb_info%bond_a, amb_info%bond_b, &
824  amb_info%bond_k, amb_info%bond_r0, particle_set, nsize, &
825  nbonh, ibh, jbh, icbh, rk, req)
826  ! Bonds non-containing hydrogens
827  CALL post_process_bonds_info(amb_info%bond_a, amb_info%bond_b, &
828  amb_info%bond_k, amb_info%bond_r0, particle_set, nsize, &
829  nbona, ib, jb, icb, rk, req)
830  ! Shrink arrays size to the minimal request
831  CALL reallocate(amb_info%bond_a, 1, nsize)
832  CALL reallocate(amb_info%bond_b, 1, nsize)
833  CALL reallocate(amb_info%bond_k, 1, nsize)
834  CALL reallocate(amb_info%bond_r0, 1, nsize)
835 
836  ! ----------------------------------------------------------
837  ! Force Fields informations related to bends
838  ! ----------------------------------------------------------
839  CALL reallocate(amb_info%bend_a, 1, buffer_size)
840  CALL reallocate(amb_info%bend_b, 1, buffer_size)
841  CALL reallocate(amb_info%bend_c, 1, buffer_size)
842  CALL reallocate(amb_info%bend_k, 1, buffer_size)
843  CALL reallocate(amb_info%bend_theta0, 1, buffer_size)
844  nsize = 0
845  ! Bends containing hydrogens
846  CALL post_process_bends_info(amb_info%bend_a, amb_info%bend_b, &
847  amb_info%bend_c, amb_info%bend_k, amb_info%bend_theta0, &
848  particle_set, nsize, ntheth, ith, jth, kth, icth, tk, teq)
849  ! Bends non-containing hydrogens
850  CALL post_process_bends_info(amb_info%bend_a, amb_info%bend_b, &
851  amb_info%bend_c, amb_info%bend_k, amb_info%bend_theta0, &
852  particle_set, nsize, ntheta, it, jt, kt, ict, tk, teq)
853  ! Shrink arrays size to the minimal request
854  CALL reallocate(amb_info%bend_a, 1, nsize)
855  CALL reallocate(amb_info%bend_b, 1, nsize)
856  CALL reallocate(amb_info%bend_c, 1, nsize)
857  CALL reallocate(amb_info%bend_k, 1, nsize)
858  CALL reallocate(amb_info%bend_theta0, 1, nsize)
859 
860  ! ----------------------------------------------------------
861  ! Force Fields informations related to torsions
862  ! in amb_info%phi0 we store PHI0
863  ! ----------------------------------------------------------
864 
865  CALL reallocate(amb_info%torsion_a, 1, buffer_size)
866  CALL reallocate(amb_info%torsion_b, 1, buffer_size)
867  CALL reallocate(amb_info%torsion_c, 1, buffer_size)
868  CALL reallocate(amb_info%torsion_d, 1, buffer_size)
869  CALL reallocate(amb_info%torsion_k, 1, buffer_size)
870  CALL reallocate(amb_info%torsion_m, 1, buffer_size)
871  CALL reallocate(amb_info%torsion_phi0, 1, buffer_size)
872  nsize = 0
873  ! Torsions containing hydrogens
874  CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b, &
875  amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, &
876  amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize, &
877  nphih, iph, jph, kph, lph, icph, pk, pn, phase)
878  ! Torsions non-containing hydrogens
879  CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b, &
880  amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, &
881  amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize, &
882  nphia, ip, jp, kp, lp, icp, pk, pn, phase)
883  ! Shrink arrays size to the minimal request
884  CALL reallocate(amb_info%torsion_a, 1, nsize)
885  CALL reallocate(amb_info%torsion_b, 1, nsize)
886  CALL reallocate(amb_info%torsion_c, 1, nsize)
887  CALL reallocate(amb_info%torsion_d, 1, nsize)
888  CALL reallocate(amb_info%torsion_k, 1, nsize)
889  CALL reallocate(amb_info%torsion_m, 1, nsize)
890  CALL reallocate(amb_info%torsion_phi0, 1, nsize)
891 
892  ! Sort dihedral metadata for faster lookup
893  IF (nphih + nphia /= 0) THEN
894  ALLOCATE (iwork(nphih + nphia))
895  CALL sort(amb_info%raw_torsion_id, 1, nphih + nphia, 1, 5, iwork)
896  DEALLOCATE (iwork)
897  END IF
898 
899  ! ----------------------------------------------------------
900  ! Post process of LJ parameters
901  ! ----------------------------------------------------------
902  CALL reallocate(amb_info%nonbond_a, 1, buffer_size)
903  CALL reallocate(amb_info%nonbond_eps, 1, buffer_size)
904  CALL reallocate(amb_info%nonbond_rmin2, 1, buffer_size)
905 
906  nsize = 0
907  CALL post_process_lj_info(amb_info%nonbond_a, amb_info%nonbond_eps, &
908  amb_info%nonbond_rmin2, particle_set, ntypes, nsize, iac, ico, &
909  cn1, cn2, natom)
910 
911  ! Shrink arrays size to the minimal request
912  CALL reallocate(amb_info%nonbond_a, 1, nsize)
913  CALL reallocate(amb_info%nonbond_eps, 1, nsize)
914  CALL reallocate(amb_info%nonbond_rmin2, 1, nsize)
915 
916  ! Deallocate at the end of the dirty job
917  DEALLOCATE (iac)
918  DEALLOCATE (ico)
919  DEALLOCATE (rk)
920  DEALLOCATE (req)
921  DEALLOCATE (tk)
922  DEALLOCATE (teq)
923  DEALLOCATE (pk)
924  DEALLOCATE (pn)
925  DEALLOCATE (phase)
926  DEALLOCATE (cn1)
927  DEALLOCATE (cn2)
928  DEALLOCATE (asol)
929  DEALLOCATE (bsol)
930  CALL timestop(handle2)
931  END IF
932  ! Always Deallocate
933  DEALLOCATE (ibh)
934  DEALLOCATE (jbh)
935  DEALLOCATE (icbh)
936  DEALLOCATE (ib)
937  DEALLOCATE (jb)
938  DEALLOCATE (icb)
939  DEALLOCATE (ith)
940  DEALLOCATE (jth)
941  DEALLOCATE (kth)
942  DEALLOCATE (icth)
943  DEALLOCATE (it)
944  DEALLOCATE (jt)
945  DEALLOCATE (kt)
946  DEALLOCATE (ict)
947  DEALLOCATE (iph)
948  DEALLOCATE (jph)
949  DEALLOCATE (kph)
950  DEALLOCATE (lph)
951  DEALLOCATE (icph)
952  DEALLOCATE (ip)
953  DEALLOCATE (jp)
954  DEALLOCATE (kp)
955  DEALLOCATE (lp)
956  DEALLOCATE (icp)
957  CALL parser_release(parser)
958  CALL timestop(handle)
959  RETURN
960  ! Output info Format
961 1000 FORMAT(t2, &
962  /' NATOM = ', i7, ' NTYPES = ', i7, ' NBONH = ', i7, ' MBONA = ', i7, &
963  /' NTHETH = ', i7, ' MTHETA = ', i7, ' NPHIH = ', i7, ' MPHIA = ', i7, &
964  /' NHPARM = ', i7, ' NPARM = ', i7, ' NNB = ', i7, ' NRES = ', i7, &
965  /' NBONA = ', i7, ' NTHETA = ', i7, ' NPHIA = ', i7, ' NUMBND = ', i7, &
966  /' NUMANG = ', i7, ' NPTRA = ', i7, ' NATYP = ', i7, ' NPHB = ', i7, &
967  /' IFBOX = ', i7, ' NMXRS = ', i7, ' IFCAP = ', i7, ' NEXTRA = ', i7,/)
968  END SUBROUTINE rdparm_amber_8
969 
970 ! **************************************************************************************************
971 !> \brief Low level routine to identify and rename unique atom types
972 !> \param isymbl ...
973 !> \param iwork ...
974 !> \param i ...
975 !> \param istart ...
976 !> \param charges ...
977 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
978 ! **************************************************************************************************
979  SUBROUTINE conform_atom_type_low(isymbl, iwork, i, istart, charges)
980  CHARACTER(LEN=default_string_length), DIMENSION(:) :: isymbl
981  INTEGER, DIMENSION(:) :: iwork
982  INTEGER, INTENT(IN) :: i
983  INTEGER, INTENT(INOUT) :: istart
984  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: charges
985 
986  CHARACTER(len=*), PARAMETER :: routinen = 'conform_atom_type_low'
987 
988  INTEGER :: counter, gind, handle, iend, ind, isize, &
989  j, k, kend, kstart
990  INTEGER, DIMENSION(:), POINTER :: cindx, lindx
991  REAL(kind=dp) :: ctmp
992  REAL(kind=dp), DIMENSION(:), POINTER :: cwork
993 
994  CALL timeset(routinen, handle)
995  iend = i - 1
996  isize = iend - istart + 1
997  ALLOCATE (cwork(isize))
998  ALLOCATE (lindx(isize))
999  ALLOCATE (cindx(isize))
1000  ind = 0
1001  DO k = istart, iend
1002  ind = ind + 1
1003  cwork(ind) = charges(iwork(k))
1004  lindx(ind) = k
1005  END DO
1006  CALL sort(cwork, isize, cindx)
1007 
1008  ctmp = cwork(1)
1009  counter = 1
1010  DO k = 2, isize
1011  IF (cwork(k) /= ctmp) THEN
1012  counter = counter + 1
1013  ctmp = cwork(k)
1014  END IF
1015  END DO
1016  IF (counter /= 1) THEN
1017  counter = 1
1018  kstart = 1
1019  ctmp = cwork(1)
1020  DO k = 2, isize
1021  IF (cwork(k) /= ctmp) THEN
1022  kend = k - 1
1023  DO j = kstart, kend
1024  gind = lindx(cindx(j))
1025  isymbl(gind) = trim(isymbl(gind))//adjustl(cp_to_string(counter))
1026  END DO
1027  counter = counter + 1
1028  ctmp = cwork(k)
1029  kstart = k
1030  END IF
1031  END DO
1032  kend = k - 1
1033  DO j = kstart, kend
1034  gind = lindx(cindx(j))
1035  isymbl(gind) = trim(isymbl(gind))//adjustl(cp_to_string(counter))
1036  END DO
1037  END IF
1038  DEALLOCATE (cwork)
1039  DEALLOCATE (lindx)
1040  DEALLOCATE (cindx)
1041  CALL timestop(handle)
1042  END SUBROUTINE conform_atom_type_low
1043 
1044 ! **************************************************************************************************
1045 !> \brief Set of Low level subroutines reading section for parmtop
1046 !> reading 1 array of integers of length dim
1047 !> \param parser ...
1048 !> \param section ...
1049 !> \param array1 ...
1050 !> \param dim ...
1051 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1052 ! **************************************************************************************************
1053  SUBROUTINE rd_amber_section_i1(parser, section, array1, dim)
1054  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1055  CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1056  INTEGER, DIMENSION(:) :: array1
1057  INTEGER, INTENT(IN) :: dim
1058 
1059  INTEGER :: i
1060  LOGICAL :: my_end
1061 
1062  CALL parser_get_next_line(parser, 1, at_end=my_end)
1063  i = 1
1064  DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1065  IF (parser_test_next_token(parser) == "EOL") &
1066  CALL parser_get_next_line(parser, 1, at_end=my_end)
1067  IF (my_end) EXIT
1068  CALL parser_get_object(parser, array1(i))
1069  i = i + 1
1070  END DO
1071  ! Trigger end of file aborting
1072  IF (my_end .AND. (i <= dim)) &
1073  CALL cp_abort(__location__, &
1074  "End of file while reading section "//trim(section)//" in amber topology file!")
1075  END SUBROUTINE rd_amber_section_i1
1076 
1077 ! **************************************************************************************************
1078 !> \brief Set of Low level subroutines reading section for parmtop
1079 !> reading 3 arrays of integers of length dim
1080 !> \param parser ...
1081 !> \param section ...
1082 !> \param array1 ...
1083 !> \param array2 ...
1084 !> \param array3 ...
1085 !> \param dim ...
1086 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1087 ! **************************************************************************************************
1088  SUBROUTINE rd_amber_section_i3(parser, section, array1, array2, array3, dim)
1089  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1090  CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1091  INTEGER, DIMENSION(:) :: array1, array2, array3
1092  INTEGER, INTENT(IN) :: dim
1093 
1094  INTEGER :: i
1095  LOGICAL :: my_end
1096 
1097  CALL parser_get_next_line(parser, 1, at_end=my_end)
1098  i = 1
1099  DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1100  !array1
1101  IF (parser_test_next_token(parser) == "EOL") &
1102  CALL parser_get_next_line(parser, 1, at_end=my_end)
1103  IF (my_end) EXIT
1104  CALL parser_get_object(parser, array1(i))
1105  !array2
1106  IF (parser_test_next_token(parser) == "EOL") &
1107  CALL parser_get_next_line(parser, 1, at_end=my_end)
1108  IF (my_end) EXIT
1109  CALL parser_get_object(parser, array2(i))
1110  !array3
1111  IF (parser_test_next_token(parser) == "EOL") &
1112  CALL parser_get_next_line(parser, 1, at_end=my_end)
1113  IF (my_end) EXIT
1114  CALL parser_get_object(parser, array3(i))
1115  i = i + 1
1116  END DO
1117  ! Trigger end of file aborting
1118  IF (my_end .AND. (i <= dim)) &
1119  CALL cp_abort(__location__, &
1120  "End of file while reading section "//trim(section)//" in amber topology file!")
1121  END SUBROUTINE rd_amber_section_i3
1122 
1123 ! **************************************************************************************************
1124 !> \brief Set of Low level subroutines reading section for parmtop
1125 !> reading 4 arrays of integers of length dim
1126 !> \param parser ...
1127 !> \param section ...
1128 !> \param array1 ...
1129 !> \param array2 ...
1130 !> \param array3 ...
1131 !> \param array4 ...
1132 !> \param dim ...
1133 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1134 ! **************************************************************************************************
1135  SUBROUTINE rd_amber_section_i4(parser, section, array1, array2, array3, array4, dim)
1136  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1137  CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1138  INTEGER, DIMENSION(:) :: array1, array2, array3, array4
1139  INTEGER, INTENT(IN) :: dim
1140 
1141  INTEGER :: i
1142  LOGICAL :: my_end
1143 
1144  CALL parser_get_next_line(parser, 1, at_end=my_end)
1145  i = 1
1146  DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1147  !array1
1148  IF (parser_test_next_token(parser) == "EOL") &
1149  CALL parser_get_next_line(parser, 1, at_end=my_end)
1150  IF (my_end) EXIT
1151  CALL parser_get_object(parser, array1(i))
1152  !array2
1153  IF (parser_test_next_token(parser) == "EOL") &
1154  CALL parser_get_next_line(parser, 1, at_end=my_end)
1155  IF (my_end) EXIT
1156  CALL parser_get_object(parser, array2(i))
1157  !array3
1158  IF (parser_test_next_token(parser) == "EOL") &
1159  CALL parser_get_next_line(parser, 1, at_end=my_end)
1160  IF (my_end) EXIT
1161  CALL parser_get_object(parser, array3(i))
1162  !array4
1163  IF (parser_test_next_token(parser) == "EOL") &
1164  CALL parser_get_next_line(parser, 1, at_end=my_end)
1165  IF (my_end) EXIT
1166  CALL parser_get_object(parser, array4(i))
1167  i = i + 1
1168  END DO
1169  ! Trigger end of file aborting
1170  IF (my_end .AND. (i <= dim)) &
1171  CALL cp_abort(__location__, &
1172  "End of file while reading section "//trim(section)//" in amber topology file!")
1173  END SUBROUTINE rd_amber_section_i4
1174 
1175 ! **************************************************************************************************
1176 !> \brief Set of Low level subroutines reading section for parmtop
1177 !> reading 5 arrays of integers of length dim
1178 !> \param parser ...
1179 !> \param section ...
1180 !> \param array1 ...
1181 !> \param array2 ...
1182 !> \param array3 ...
1183 !> \param array4 ...
1184 !> \param array5 ...
1185 !> \param dim ...
1186 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1187 ! **************************************************************************************************
1188  SUBROUTINE rd_amber_section_i5(parser, section, array1, array2, array3, array4, &
1189  array5, dim)
1190  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1191  CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1192  INTEGER, DIMENSION(:) :: array1, array2, array3, array4, array5
1193  INTEGER, INTENT(IN) :: dim
1194 
1195  INTEGER :: i
1196  LOGICAL :: my_end
1197 
1198  CALL parser_get_next_line(parser, 1, at_end=my_end)
1199  i = 1
1200  DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1201  !array1
1202  IF (parser_test_next_token(parser) == "EOL") &
1203  CALL parser_get_next_line(parser, 1, at_end=my_end)
1204  IF (my_end) EXIT
1205  CALL parser_get_object(parser, array1(i))
1206  !array2
1207  IF (parser_test_next_token(parser) == "EOL") &
1208  CALL parser_get_next_line(parser, 1, at_end=my_end)
1209  IF (my_end) EXIT
1210  CALL parser_get_object(parser, array2(i))
1211  !array3
1212  IF (parser_test_next_token(parser) == "EOL") &
1213  CALL parser_get_next_line(parser, 1, at_end=my_end)
1214  IF (my_end) EXIT
1215  CALL parser_get_object(parser, array3(i))
1216  !array4
1217  IF (parser_test_next_token(parser) == "EOL") &
1218  CALL parser_get_next_line(parser, 1, at_end=my_end)
1219  IF (my_end) EXIT
1220  CALL parser_get_object(parser, array4(i))
1221  !array5
1222  IF (parser_test_next_token(parser) == "EOL") &
1223  CALL parser_get_next_line(parser, 1, at_end=my_end)
1224  IF (my_end) EXIT
1225  CALL parser_get_object(parser, array5(i))
1226  i = i + 1
1227  END DO
1228  ! Trigger end of file aborting
1229  IF (my_end .AND. (i <= dim)) &
1230  CALL cp_abort(__location__, &
1231  "End of file while reading section "//trim(section)//" in amber topology file!")
1232  END SUBROUTINE rd_amber_section_i5
1233 
1234 ! **************************************************************************************************
1235 !> \brief Set of Low level subroutines reading section for parmtop
1236 !> reading 1 array of strings of length dim
1237 !> \param parser ...
1238 !> \param section ...
1239 !> \param array1 ...
1240 !> \param dim ...
1241 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1242 ! **************************************************************************************************
1243  SUBROUTINE rd_amber_section_c1(parser, section, array1, dim)
1244  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1245  CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1246  CHARACTER(LEN=default_string_length), DIMENSION(:) :: array1
1247  INTEGER, INTENT(IN) :: dim
1248 
1249  INTEGER :: i
1250  LOGICAL :: my_end
1251 
1252  CALL parser_get_next_line(parser, 1, at_end=my_end)
1253  i = 1
1254  DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1255  IF (parser_test_next_token(parser) == "EOL") &
1256  CALL parser_get_next_line(parser, 1, at_end=my_end)
1257  IF (my_end) EXIT
1258  CALL parser_get_object(parser, array1(i), lower_to_upper=.true.)
1259  i = i + 1
1260  END DO
1261  ! Trigger end of file aborting
1262  IF (my_end .AND. (i <= dim)) &
1263  CALL cp_abort(__location__, &
1264  "End of file while reading section "//trim(section)//" in amber topology file!")
1265  END SUBROUTINE rd_amber_section_c1
1266 
1267 ! **************************************************************************************************
1268 !> \brief Set of Low level subroutines reading section for parmtop
1269 !> reading 1 array of strings of length dim
1270 !> \param parser ...
1271 !> \param section ...
1272 !> \param array1 ...
1273 !> \param dim ...
1274 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1275 ! **************************************************************************************************
1276  SUBROUTINE rd_amber_section_r1(parser, section, array1, dim)
1277  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1278  CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1279  REAL(kind=dp), DIMENSION(:) :: array1
1280  INTEGER, INTENT(IN) :: dim
1281 
1282  INTEGER :: i
1283  LOGICAL :: my_end
1284 
1285  CALL parser_get_next_line(parser, 1, at_end=my_end)
1286  i = 1
1287  DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1288  IF (parser_test_next_token(parser) == "EOL") &
1289  CALL parser_get_next_line(parser, 1, at_end=my_end)
1290  IF (my_end) EXIT
1291  CALL parser_get_object(parser, array1(i))
1292  i = i + 1
1293  END DO
1294  ! Trigger end of file aborting
1295  IF (my_end .AND. (i <= dim)) &
1296  CALL cp_abort(__location__, &
1297  "End of file while reading section "//trim(section)//" in amber topology file!")
1298  END SUBROUTINE rd_amber_section_r1
1299 
1300 ! **************************************************************************************************
1301 !> \brief Check the version of the AMBER topology file (we can handle from v8 on)
1302 !> \param parser ...
1303 !> \param section ...
1304 !> \param input_format ...
1305 !> \return ...
1306 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1307 ! **************************************************************************************************
1308  FUNCTION get_section_parmtop(parser, section, input_format) RESULT(another_section)
1309  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1310  CHARACTER(LEN=default_string_length), INTENT(OUT) :: section, input_format
1311  LOGICAL :: another_section
1312 
1313  INTEGER :: end_f, indflag, start_f
1314  LOGICAL :: found, my_end
1315 
1316  CALL parser_search_string(parser, "%FLAG", .true., found, begin_line=.true.)
1317  IF (found) THEN
1318  ! section label
1319  indflag = index(parser%input_line, "%FLAG") + len_trim("%FLAG")
1320  DO WHILE (index(parser%input_line(indflag:indflag), " ") /= 0)
1321  indflag = indflag + 1
1322  END DO
1323  section = trim(parser%input_line(indflag:))
1324  ! Input format
1325  CALL parser_get_next_line(parser, 1, at_end=my_end)
1326  IF (index(parser%input_line, "%FORMAT") == 0 .OR. my_end) &
1327  cpabort("Expecting %FORMAT. Not found! Abort reading of AMBER topology file!")
1328 
1329  start_f = index(parser%input_line, "(")
1330  end_f = index(parser%input_line, ")")
1331  input_format = parser%input_line(start_f:end_f)
1332  another_section = .true.
1333  ELSE
1334  another_section = .false.
1335  END IF
1336  END FUNCTION get_section_parmtop
1337 
1338 ! **************************************************************************************************
1339 !> \brief Check the version of the AMBER topology file (we can handle from v8 on)
1340 !> \param parser ...
1341 !> \param output_unit ...
1342 !> \return ...
1343 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1344 ! **************************************************************************************************
1345  FUNCTION check_amber_8_std(parser, output_unit) RESULT(found_AMBER_V8)
1346  TYPE(cp_parser_type), INTENT(INOUT) :: parser
1347  INTEGER, INTENT(IN) :: output_unit
1348  LOGICAL :: found_amber_v8
1349 
1350  CALL parser_search_string(parser, "%VERSION ", .true., found_amber_v8, begin_line=.true.)
1351  IF (.NOT. found_amber_v8) &
1352  CALL cp_abort(__location__, &
1353  "This is not an AMBER V.8 PRMTOP format file. Cannot interpret older "// &
1354  "AMBER file formats. ")
1355  IF (output_unit > 0) WRITE (output_unit, '(" AMBER_INFO| ",A)') "Amber PrmTop V.8 or greater.", &
1356  trim(parser%input_line)
1357 
1358  END FUNCTION check_amber_8_std
1359 
1360 ! **************************************************************************************************
1361 !> \brief Post processing of forcefield information related to bonds
1362 !> \param label_a ...
1363 !> \param label_b ...
1364 !> \param k ...
1365 !> \param r0 ...
1366 !> \param particle_set ...
1367 !> \param ibond ...
1368 !> \param nbond ...
1369 !> \param ib ...
1370 !> \param jb ...
1371 !> \param icb ...
1372 !> \param rk ...
1373 !> \param req ...
1374 !> \author Teodoro Laino [tlaino] - 11.2008
1375 ! **************************************************************************************************
1376  SUBROUTINE post_process_bonds_info(label_a, label_b, k, r0, particle_set, ibond, &
1377  nbond, ib, jb, icb, rk, req)
1378  CHARACTER(LEN=default_string_length), &
1379  DIMENSION(:), POINTER :: label_a, label_b
1380  REAL(kind=dp), DIMENSION(:), POINTER :: k, r0
1381  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1382  INTEGER, INTENT(INOUT) :: ibond
1383  INTEGER, INTENT(IN) :: nbond
1384  INTEGER, DIMENSION(:), INTENT(IN) :: ib, jb, icb
1385  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rk, req
1386 
1387  CHARACTER(len=*), PARAMETER :: routinen = 'post_process_bonds_info'
1388 
1389  CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b
1390  CHARACTER(LEN=default_string_length), &
1391  ALLOCATABLE, DIMENSION(:, :) :: work_label
1392  INTEGER :: handle, i
1393  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1394  LOGICAL :: l_dum
1395 
1396  CALL timeset(routinen, handle)
1397  IF (nbond /= 0) THEN
1398  ALLOCATE (work_label(2, nbond))
1399  ALLOCATE (iwork(nbond))
1400  DO i = 1, nbond
1401  name_atm_a = particle_set(ib(i))%atomic_kind%name
1402  name_atm_b = particle_set(jb(i))%atomic_kind%name
1403  l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
1404  work_label(1, i) = name_atm_a
1405  work_label(2, i) = name_atm_b
1406  END DO
1407  CALL sort(work_label, 1, nbond, 1, 2, iwork)
1408 
1409  ibond = ibond + 1
1410  ! In case we need more space ... give it up...
1411  IF (ibond > SIZE(label_a)) THEN
1412  CALL reallocate(label_a, 1, int(buffer_size + ibond*1.5_dp))
1413  CALL reallocate(label_b, 1, int(buffer_size + ibond*1.5_dp))
1414  CALL reallocate(k, 1, int(buffer_size + ibond*1.5_dp))
1415  CALL reallocate(r0, 1, int(buffer_size + ibond*1.5_dp))
1416  END IF
1417  label_a(ibond) = work_label(1, 1)
1418  label_b(ibond) = work_label(2, 1)
1419  k(ibond) = rk(icb(iwork(1)))
1420  r0(ibond) = req(icb(iwork(1)))
1421 
1422  DO i = 2, nbond
1423  IF ((work_label(1, i) /= label_a(ibond)) .OR. &
1424  (work_label(2, i) /= label_b(ibond))) THEN
1425  ibond = ibond + 1
1426  ! In case we need more space ... give it up...
1427  IF (ibond > SIZE(label_a)) THEN
1428  CALL reallocate(label_a, 1, int(buffer_size + ibond*1.5_dp))
1429  CALL reallocate(label_b, 1, int(buffer_size + ibond*1.5_dp))
1430  CALL reallocate(k, 1, int(buffer_size + ibond*1.5_dp))
1431  CALL reallocate(r0, 1, int(buffer_size + ibond*1.5_dp))
1432  END IF
1433  label_a(ibond) = work_label(1, i)
1434  label_b(ibond) = work_label(2, i)
1435  k(ibond) = rk(icb(iwork(i)))
1436  r0(ibond) = req(icb(iwork(i)))
1437  END IF
1438  END DO
1439 
1440  DEALLOCATE (work_label)
1441  DEALLOCATE (iwork)
1442  END IF
1443  CALL timestop(handle)
1444  END SUBROUTINE post_process_bonds_info
1445 
1446 ! **************************************************************************************************
1447 !> \brief Post processing of forcefield information related to bends
1448 !> \param label_a ...
1449 !> \param label_b ...
1450 !> \param label_c ...
1451 !> \param k ...
1452 !> \param theta0 ...
1453 !> \param particle_set ...
1454 !> \param itheta ...
1455 !> \param ntheta ...
1456 !> \param it ...
1457 !> \param jt ...
1458 !> \param kt ...
1459 !> \param ict ...
1460 !> \param tk ...
1461 !> \param teq ...
1462 !> \author Teodoro Laino [tlaino] - 11.2008
1463 ! **************************************************************************************************
1464  SUBROUTINE post_process_bends_info(label_a, label_b, label_c, k, theta0, &
1465  particle_set, itheta, ntheta, it, jt, kt, ict, tk, teq)
1466  CHARACTER(LEN=default_string_length), &
1467  DIMENSION(:), POINTER :: label_a, label_b, label_c
1468  REAL(kind=dp), DIMENSION(:), POINTER :: k, theta0
1469  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1470  INTEGER, INTENT(INOUT) :: itheta
1471  INTEGER, INTENT(IN) :: ntheta
1472  INTEGER, DIMENSION(:), INTENT(IN) :: it, jt, kt, ict
1473  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: tk, teq
1474 
1475  CHARACTER(len=*), PARAMETER :: routinen = 'post_process_bends_info'
1476 
1477  CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c
1478  CHARACTER(LEN=default_string_length), &
1479  ALLOCATABLE, DIMENSION(:, :) :: work_label
1480  INTEGER :: handle, i
1481  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1482  LOGICAL :: l_dum
1483 
1484  CALL timeset(routinen, handle)
1485  IF (ntheta /= 0) THEN
1486  ALLOCATE (work_label(3, ntheta))
1487  ALLOCATE (iwork(ntheta))
1488  DO i = 1, ntheta
1489  name_atm_a = particle_set(it(i))%atomic_kind%name
1490  name_atm_b = particle_set(jt(i))%atomic_kind%name
1491  name_atm_c = particle_set(kt(i))%atomic_kind%name
1492  l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, &
1493  id3=name_atm_c)
1494  work_label(1, i) = name_atm_a
1495  work_label(2, i) = name_atm_b
1496  work_label(3, i) = name_atm_c
1497  END DO
1498 
1499  CALL sort(work_label, 1, ntheta, 1, 3, iwork)
1500 
1501  itheta = itheta + 1
1502  ! In case we need more space ... give it up...
1503  IF (itheta > SIZE(label_a)) THEN
1504  CALL reallocate(label_a, 1, int(buffer_size + itheta*1.5_dp))
1505  CALL reallocate(label_b, 1, int(buffer_size + itheta*1.5_dp))
1506  CALL reallocate(label_c, 1, int(buffer_size + itheta*1.5_dp))
1507  CALL reallocate(k, 1, int(buffer_size + itheta*1.5_dp))
1508  CALL reallocate(theta0, 1, int(buffer_size + itheta*1.5_dp))
1509  END IF
1510  label_a(itheta) = work_label(1, 1)
1511  label_b(itheta) = work_label(2, 1)
1512  label_c(itheta) = work_label(3, 1)
1513  k(itheta) = tk(ict(iwork(1)))
1514  theta0(itheta) = teq(ict(iwork(1)))
1515 
1516  DO i = 2, ntheta
1517  IF ((work_label(1, i) /= label_a(itheta)) .OR. &
1518  (work_label(2, i) /= label_b(itheta)) .OR. &
1519  (work_label(3, i) /= label_c(itheta))) THEN
1520  itheta = itheta + 1
1521  ! In case we need more space ... give it up...
1522  IF (itheta > SIZE(label_a)) THEN
1523  CALL reallocate(label_a, 1, int(buffer_size + itheta*1.5_dp))
1524  CALL reallocate(label_b, 1, int(buffer_size + itheta*1.5_dp))
1525  CALL reallocate(label_c, 1, int(buffer_size + itheta*1.5_dp))
1526  CALL reallocate(k, 1, int(buffer_size + itheta*1.5_dp))
1527  CALL reallocate(theta0, 1, int(buffer_size + itheta*1.5_dp))
1528  END IF
1529  label_a(itheta) = work_label(1, i)
1530  label_b(itheta) = work_label(2, i)
1531  label_c(itheta) = work_label(3, i)
1532  k(itheta) = tk(ict(iwork(i)))
1533  theta0(itheta) = teq(ict(iwork(i)))
1534  END IF
1535  END DO
1536 
1537  DEALLOCATE (work_label)
1538  DEALLOCATE (iwork)
1539  END IF
1540  CALL timestop(handle)
1541  END SUBROUTINE post_process_bends_info
1542 
1543 ! **************************************************************************************************
1544 !> \brief Post processing of forcefield information related to torsions
1545 !> \param label_a ...
1546 !> \param label_b ...
1547 !> \param label_c ...
1548 !> \param label_d ...
1549 !> \param k ...
1550 !> \param m ...
1551 !> \param phi0 ...
1552 !> \param particle_set ...
1553 !> \param iphi ...
1554 !> \param nphi ...
1555 !> \param ip ...
1556 !> \param jp ...
1557 !> \param kp ...
1558 !> \param lp ...
1559 !> \param icp ...
1560 !> \param pk ...
1561 !> \param pn ...
1562 !> \param phase ...
1563 !> \author Teodoro Laino [tlaino] - 11.2008
1564 ! **************************************************************************************************
1565  SUBROUTINE post_process_torsions_info(label_a, label_b, label_c, label_d, k, &
1566  m, phi0, particle_set, iphi, nphi, ip, jp, kp, lp, icp, pk, pn, phase)
1567  CHARACTER(LEN=default_string_length), &
1568  DIMENSION(:), POINTER :: label_a, label_b, label_c, label_d
1569  REAL(kind=dp), DIMENSION(:), POINTER :: k
1570  INTEGER, DIMENSION(:), POINTER :: m
1571  REAL(kind=dp), DIMENSION(:), POINTER :: phi0
1572  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1573  INTEGER, INTENT(INOUT) :: iphi
1574  INTEGER, INTENT(IN) :: nphi
1575  INTEGER, DIMENSION(:), INTENT(IN) :: ip, jp, kp, lp, icp
1576  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: pk, pn, phase
1577 
1578  CHARACTER(len=*), PARAMETER :: routinen = 'post_process_torsions_info'
1579 
1580  CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c, &
1581  name_atm_d
1582  CHARACTER(LEN=default_string_length), &
1583  ALLOCATABLE, DIMENSION(:, :) :: work_label
1584  INTEGER :: handle, i
1585  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1586  LOGICAL :: l_dum
1587 
1588  CALL timeset(routinen, handle)
1589  IF (nphi /= 0) THEN
1590  ALLOCATE (work_label(6, nphi))
1591  ALLOCATE (iwork(nphi))
1592  DO i = 1, nphi
1593  name_atm_a = particle_set(ip(i))%atomic_kind%name
1594  name_atm_b = particle_set(jp(i))%atomic_kind%name
1595  name_atm_c = particle_set(kp(i))%atomic_kind%name
1596  name_atm_d = particle_set(lp(i))%atomic_kind%name
1597  l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, &
1598  id3=name_atm_c, id4=name_atm_d)
1599  work_label(1, i) = name_atm_a
1600  work_label(2, i) = name_atm_b
1601  work_label(3, i) = name_atm_c
1602  work_label(4, i) = name_atm_d
1603  ! Phase and multiplicity must be kept into account
1604  ! for the ordering of the torsions
1605  work_label(5, i) = trim(adjustl(cp_to_string(phase(icp(i)))))
1606  work_label(6, i) = trim(adjustl(cp_to_string(pn(icp(i)))))
1607  END DO
1608 
1609  CALL sort(work_label, 1, nphi, 1, 6, iwork)
1610 
1611  iphi = iphi + 1
1612  ! In case we need more space ... give it up...
1613  IF (iphi > SIZE(label_a)) THEN
1614  CALL reallocate(label_a, 1, int(buffer_size + iphi*1.5_dp))
1615  CALL reallocate(label_b, 1, int(buffer_size + iphi*1.5_dp))
1616  CALL reallocate(label_c, 1, int(buffer_size + iphi*1.5_dp))
1617  CALL reallocate(label_d, 1, int(buffer_size + iphi*1.5_dp))
1618  CALL reallocate(k, 1, int(buffer_size + iphi*1.5_dp))
1619  CALL reallocate(m, 1, int(buffer_size + iphi*1.5_dp))
1620  CALL reallocate(phi0, 1, int(buffer_size + iphi*1.5_dp))
1621  END IF
1622  label_a(iphi) = work_label(1, 1)
1623  label_b(iphi) = work_label(2, 1)
1624  label_c(iphi) = work_label(3, 1)
1625  label_d(iphi) = work_label(4, 1)
1626  k(iphi) = pk(icp(iwork(1)))
1627  m(iphi) = nint(pn(icp(iwork(1))))
1628  IF (m(iphi) - pn(icp(iwork(1))) .GT. epsilon(1.0_dp)) THEN
1629  ! non integer torsions not supported
1630  cpabort("")
1631  END IF
1632 
1633  phi0(iphi) = phase(icp(iwork(1)))
1634 
1635  DO i = 2, nphi
1636  ! We don't consider the possibility that a torsion can have same
1637  ! phase, periodicity but different value of k.. in this case the
1638  ! potential should be summed-up
1639  IF ((work_label(1, i) /= label_a(iphi)) .OR. &
1640  (work_label(2, i) /= label_b(iphi)) .OR. &
1641  (work_label(3, i) /= label_c(iphi)) .OR. &
1642  (work_label(4, i) /= label_d(iphi)) .OR. &
1643  (pn(icp(iwork(i))) /= m(iphi)) .OR. &
1644  (phase(icp(iwork(i))) /= phi0(iphi))) THEN
1645  iphi = iphi + 1
1646  ! In case we need more space ... give it up...
1647  IF (iphi > SIZE(label_a)) THEN
1648  CALL reallocate(label_a, 1, int(buffer_size + iphi*1.5_dp))
1649  CALL reallocate(label_b, 1, int(buffer_size + iphi*1.5_dp))
1650  CALL reallocate(label_c, 1, int(buffer_size + iphi*1.5_dp))
1651  CALL reallocate(label_d, 1, int(buffer_size + iphi*1.5_dp))
1652  CALL reallocate(k, 1, int(buffer_size + iphi*1.5_dp))
1653  CALL reallocate(m, 1, int(buffer_size + iphi*1.5_dp))
1654  CALL reallocate(phi0, 1, int(buffer_size + iphi*1.5_dp))
1655  END IF
1656  label_a(iphi) = work_label(1, i)
1657  label_b(iphi) = work_label(2, i)
1658  label_c(iphi) = work_label(3, i)
1659  label_d(iphi) = work_label(4, i)
1660  k(iphi) = pk(icp(iwork(i)))
1661  m(iphi) = nint(pn(icp(iwork(i))))
1662  IF (m(iphi) - pn(icp(iwork(i))) .GT. epsilon(1.0_dp)) THEN
1663  ! non integer torsions not supported
1664  cpabort("")
1665  END IF
1666  phi0(iphi) = phase(icp(iwork(i)))
1667  END IF
1668  END DO
1669 
1670  DEALLOCATE (work_label)
1671  DEALLOCATE (iwork)
1672  END IF
1673  CALL timestop(handle)
1674  END SUBROUTINE post_process_torsions_info
1675 
1676 ! **************************************************************************************************
1677 !> \brief Post processing of forcefield information related to Lennard-Jones
1678 !> \param atom_label ...
1679 !> \param eps ...
1680 !> \param sigma ...
1681 !> \param particle_set ...
1682 !> \param ntypes ...
1683 !> \param nsize ...
1684 !> \param iac ...
1685 !> \param ico ...
1686 !> \param cn1 ...
1687 !> \param cn2 ...
1688 !> \param natom ...
1689 !> \author Teodoro Laino [tlaino] - 11.2008
1690 ! **************************************************************************************************
1691  SUBROUTINE post_process_lj_info(atom_label, eps, sigma, particle_set, &
1692  ntypes, nsize, iac, ico, cn1, cn2, natom)
1693  CHARACTER(LEN=default_string_length), &
1694  DIMENSION(:), POINTER :: atom_label
1695  REAL(kind=dp), DIMENSION(:), POINTER :: eps, sigma
1696  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1697  INTEGER, INTENT(IN) :: ntypes
1698  INTEGER, INTENT(INOUT) :: nsize
1699  INTEGER, DIMENSION(:), INTENT(IN) :: iac, ico
1700  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: cn1, cn2
1701  INTEGER, INTENT(IN) :: natom
1702 
1703  CHARACTER(len=*), PARAMETER :: routinen = 'post_process_LJ_info'
1704 
1705  CHARACTER(LEN=default_string_length) :: name_atm_a
1706  CHARACTER(LEN=default_string_length), &
1707  ALLOCATABLE, DIMENSION(:) :: work_label
1708  INTEGER :: handle, i
1709  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1710  LOGICAL :: check, l_dum
1711  REAL(kind=dp) :: f12, f6, my_eps, my_sigma, sigma6
1712 
1713  CALL timeset(routinen, handle)
1714  ALLOCATE (work_label(natom))
1715  ALLOCATE (iwork(natom))
1716  DO i = 1, natom
1717  name_atm_a = particle_set(i)%atomic_kind%name
1718  l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a)
1719  work_label(i) = name_atm_a
1720  END DO
1721  CALL sort(work_label, natom, iwork)
1722 
1723  nsize = nsize + 1
1724  IF (nsize > SIZE(atom_label)) THEN
1725  CALL reallocate(atom_label, 1, int(buffer_size + nsize*1.5_dp))
1726  CALL reallocate(eps, 1, int(buffer_size + nsize*1.5_dp))
1727  CALL reallocate(sigma, 1, int(buffer_size + nsize*1.5_dp))
1728  END IF
1729  f12 = cn1(ico(ntypes*(iac(iwork(1)) - 1) + iac(iwork(1))))
1730  f6 = cn2(ico(ntypes*(iac(iwork(1)) - 1) + iac(iwork(1))))
1731  check = (f6 == 0.0_dp) .EQV. (f12 == 0.0_dp)
1732  cpassert(check)
1733  my_sigma = 0.0_dp
1734  my_eps = 0.0_dp
1735  IF (f6 /= 0.0_dp) THEN
1736  sigma6 = (2.0_dp*f12/f6)
1737  my_sigma = sigma6**(1.0_dp/6.0_dp)
1738  my_eps = f6/(2.0_dp*sigma6)
1739  END IF
1740  atom_label(nsize) = work_label(1)
1741  sigma(nsize) = my_sigma/2.0_dp
1742  eps(nsize) = my_eps
1743 
1744  DO i = 2, natom
1745  IF (work_label(i) /= atom_label(nsize)) THEN
1746  nsize = nsize + 1
1747  ! In case we need more space ... give it up...
1748  IF (nsize > SIZE(atom_label)) THEN
1749  CALL reallocate(atom_label, 1, int(buffer_size + nsize*1.5_dp))
1750  CALL reallocate(eps, 1, int(buffer_size + nsize*1.5_dp))
1751  CALL reallocate(sigma, 1, int(buffer_size + nsize*1.5_dp))
1752  END IF
1753  f12 = cn1(ico(ntypes*(iac(iwork(i)) - 1) + iac(iwork(i))))
1754  f6 = cn2(ico(ntypes*(iac(iwork(i)) - 1) + iac(iwork(i))))
1755  check = (f6 == 0.0_dp) .EQV. (f12 == 0.0_dp)
1756  cpassert(check)
1757  my_sigma = 0.0_dp
1758  my_eps = 0.0_dp
1759  IF (f6 /= 0.0_dp) THEN
1760  sigma6 = (2.0_dp*f12/f6)
1761  my_sigma = sigma6**(1.0_dp/6.0_dp)
1762  my_eps = f6/(2.0_dp*sigma6)
1763  END IF
1764  atom_label(nsize) = work_label(i)
1765  sigma(nsize) = my_sigma/2.0_dp
1766  eps(nsize) = my_eps
1767  END IF
1768  END DO
1769 
1770  DEALLOCATE (work_label)
1771  DEALLOCATE (iwork)
1772  CALL timestop(handle)
1773  END SUBROUTINE post_process_lj_info
1774 
1775 END MODULE topology_amber
1776 
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
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Definition: cp_parser_methods.F:421
cp_parser_methods::parser_test_next_token
character(len=3) function, public parser_test_next_token(parser, string_length)
Test next input object.
Definition: cp_parser_methods.F:691
cp_parser_methods::parser_search_string
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Definition: cp_parser_methods.F:786
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition: cp_parser_types.F:109
cp_parser_types::parser_create
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Definition: cp_parser_types.F:147
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition: cp_units.F:1150
force_field_types
Define all structures types related to force_fields.
Definition: force_field_types.F:16
input_cp2k_restarts_util
Definition: input_cp2k_restarts_util.F:8
input_cp2k_restarts_util::section_velocity_val_set
subroutine, public section_velocity_val_set(velocity_section, particles, velocity, conv_factor)
routine to dump velocities.. fast implementation
Definition: input_cp2k_restarts_util.F:44
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_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
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
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
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
string_table
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition: string_table.F:22
string_table::s2s
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition: string_table.F:141
string_table::str2id
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
Definition: string_table.F:72
string_table::id2str
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition: string_table.F:115
topology_amber
Handles all functions used to read and interpret AMBER coordinates and topology files.
Definition: topology_amber.F:14
topology_amber::rdparm_amber_8
subroutine, public rdparm_amber_8(filename, output_unit, para_env, do_connectivity, do_forcefield, atom_info, conn_info, amb_info, particle_set)
Access information form the AMBER topology file Notes on file structure:
Definition: topology_amber.F:368
topology_amber::read_coordinate_crd
subroutine, public read_coordinate_crd(topology, para_env, subsys_section)
Reads the ‘coord’ version generated by the PARM or LEaP programs, as well as the ‘restrt’ version,...
Definition: topology_amber.F:101
topology_amber::read_connectivity_amber
subroutine, public read_connectivity_amber(filename, topology, para_env, subsys_section)
Read AMBER topology file (.top) : At this level we parse only the connectivity info the ....
Definition: topology_amber.F:286
topology_generate_util
Collection of subroutine needed for topology related things.
Definition: topology_generate_util.F:13
topology_generate_util::topology_generate_molname
subroutine, public topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, id_molname)
Generates molnames: useful when the connectivity on file does not provide them.
Definition: topology_generate_util.F:95
topology_types
Definition: topology_types.F:18
topology
Control for reading in different topologies and coordinates.
Definition: topology.F:13
util
All kind of helpful little routines.
Definition: util.F:14
topology_types::atom_info_type
Definition: topology_types.F:35
topology_types::connectivity_info_type
Definition: topology_types.F:53
topology_types::topology_parameters_type
Definition: topology_types.F:142