(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
14MODULE 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
64CONTAINS
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
9611000 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
1775END MODULE topology_amber
1776
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_parser_methods::parser_get_object
Definition cp_parser_methods.F:54
memory_utilities::reallocate
Definition memory_utilities.F:27
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
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_parser_types::cp_parser_type
represent a parser
Definition cp_parser_types.F:77
force_field_types::amber_info_type
Definition force_field_types.F:136
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
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