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