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