Line data Source code
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 : ! **************************************************************************************************
14 : MODULE topology_amber
15 : USE cp_log_handling, ONLY: cp_get_default_logger,&
16 : cp_logger_type,&
17 : cp_to_string
18 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 : cp_print_key_unit_nr
20 : USE cp_parser_methods, ONLY: parser_get_next_line,&
21 : parser_get_object,&
22 : parser_search_string,&
23 : parser_test_next_token
24 : USE cp_parser_types, ONLY: cp_parser_type,&
25 : parser_create,&
26 : parser_release
27 : USE cp_units, ONLY: cp_unit_to_cp2k
28 : USE force_field_types, ONLY: amber_info_type
29 : USE input_cp2k_restarts_util, ONLY: section_velocity_val_set
30 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 : section_vals_type
32 : USE kinds, ONLY: default_string_length,&
33 : dp
34 : USE memory_utilities, ONLY: reallocate
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE particle_types, ONLY: particle_type
37 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
38 : USE string_table, ONLY: id2str,&
39 : s2s,&
40 : str2id
41 : USE topology_generate_util, ONLY: topology_generate_molname
42 : USE topology_types, ONLY: atom_info_type,&
43 : connectivity_info_type,&
44 : topology_parameters_type
45 : USE util, ONLY: sort
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_amber'
51 : REAL(KIND=dp), PARAMETER, PRIVATE :: amber_conv_factor = 20.4550_dp, &
52 : amber_conv_charge = 18.2223_dp
53 : INTEGER, PARAMETER, PRIVATE :: buffer_size = 1
54 :
55 : PRIVATE
56 : PUBLIC :: read_coordinate_crd, read_connectivity_amber, rdparm_amber_8
57 :
58 : ! Reading Amber sections routines
59 : INTERFACE rd_amber_section
60 : MODULE PROCEDURE rd_amber_section_i1, rd_amber_section_c1, rd_amber_section_r1, &
61 : rd_amber_section_i3, rd_amber_section_i4, rd_amber_section_i5
62 : END INTERFACE
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief Reads the `coord' version generated by the PARM or LEaP programs, as
68 : !> well as the `restrt' version, resulting from energy minimization or
69 : !> molecular dynamics in SANDER or GIBBS. It may contain velocity and
70 : !> periodic box information.
71 : !>
72 : !> Official Format from the AMBER homepage
73 : !> FORMAT(20A4) ITITL
74 : !> ITITL : the title of the current run, from the AMBER
75 : !> parameter/topology file
76 : !>
77 : !> FORMAT(I5,5E15.7) NATOM,TIME
78 : !> NATOM : total number of atoms in coordinate file
79 : !> TIME : option, current time in the simulation (picoseconds)
80 : !>
81 : !> FORMAT(6F12.7) (X(i), Y(i), Z(i), i = 1,NATOM)
82 : !> X,Y,Z : coordinates
83 : !>
84 : !> IF dynamics
85 : !>
86 : !> FORMAT(6F12.7) (VX(i), VY(i), VZ(i), i = 1,NATOM)
87 : !> VX,VY,VZ : velocities (units: Angstroms per 1/20.455 ps)
88 : !>
89 : !> IF constant pressure (in 4.1, also constant volume)
90 : !>
91 : !> FORMAT(6F12.7) BOX(1), BOX(2), BOX(3)
92 : !> BOX : size of the periodic box
93 : !>
94 : !>
95 : !> \param topology ...
96 : !> \param para_env ...
97 : !> \param subsys_section ...
98 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
99 : ! **************************************************************************************************
100 104 : 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 26 : 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 26 : NULLIFY (logger, velocity)
117 52 : logger => cp_get_default_logger()
118 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/CRD_INFO", &
119 26 : extension=".subsysLog")
120 26 : CALL timeset(routineN, handle)
121 :
122 26 : atom_info => topology%atom_info
123 26 : IF (iw > 0) WRITE (iw, *) " Reading in CRD file ", TRIM(topology%coord_file_name)
124 :
125 : ! Title Section
126 26 : IF (iw > 0) WRITE (iw, '(T2,A)') 'CRD_INFO| Parsing the TITLE section'
127 26 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
128 26 : CALL parser_get_next_line(parser, 1)
129 : ! Title may be missing
130 26 : IF (parser_test_next_token(parser) == "STR") THEN
131 20 : CALL parser_get_object(parser, string, string_length=default_string_length)
132 20 : IF (iw > 0) WRITE (iw, '(T2,A)') 'CRD_INFO| '//TRIM(string)
133 : ! Natom and Time (which we ignore)
134 46 : CALL parser_get_next_line(parser, 1)
135 : END IF
136 26 : CALL parser_get_object(parser, natom)
137 26 : topology%natoms = natom
138 26 : IF (iw > 0) WRITE (iw, '(T2,A,I0)') 'CRD_INFO| Number of atoms: ', natom
139 26 : CALL reallocate(atom_info%id_molname, 1, natom)
140 26 : CALL reallocate(atom_info%id_resname, 1, natom)
141 26 : CALL reallocate(atom_info%resid, 1, natom)
142 26 : CALL reallocate(atom_info%id_atmname, 1, natom)
143 26 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
144 26 : CALL reallocate(atom_info%atm_mass, 1, natom)
145 26 : CALL reallocate(atom_info%atm_charge, 1, natom)
146 26 : CALL reallocate(atom_info%occup, 1, natom)
147 26 : CALL reallocate(atom_info%beta, 1, natom)
148 26 : CALL reallocate(atom_info%id_element, 1, natom)
149 :
150 : ! Element is assigned on the basis of the atm_name
151 26 : topology%aa_element = .TRUE.
152 :
153 : ! Coordinates
154 26 : CALL parser_get_next_line(parser, 1, at_end=my_end)
155 38826 : DO j = 1, natom - MOD(natom, 2), 2
156 38800 : IF (my_end) EXIT
157 38800 : READ (parser%input_line, *) atom_info%r(1, j), atom_info%r(2, j), atom_info%r(3, j), &
158 77600 : 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 38800 : atom_info%id_atmname(j) = str2id(s2s("__UNDEF__"))
162 38800 : atom_info%id_molname(j) = str2id(s2s("__UNDEF__"))
163 38800 : atom_info%id_resname(j) = str2id(s2s("__UNDEF__"))
164 38800 : atom_info%id_element(j) = str2id(s2s("__UNDEF__"))
165 38800 : atom_info%resid(j) = HUGE(0)
166 38800 : atom_info%atm_mass(j) = HUGE(0.0_dp)
167 38800 : atom_info%atm_charge(j) = -HUGE(0.0_dp)
168 38800 : atom_info%r(1, j) = cp_unit_to_cp2k(atom_info%r(1, j), "angstrom")
169 38800 : atom_info%r(2, j) = cp_unit_to_cp2k(atom_info%r(2, j), "angstrom")
170 38800 : atom_info%r(3, j) = cp_unit_to_cp2k(atom_info%r(3, j), "angstrom")
171 :
172 38800 : atom_info%id_atmname(j + 1) = str2id(s2s("__UNDEF__"))
173 38800 : atom_info%id_molname(j + 1) = str2id(s2s("__UNDEF__"))
174 38800 : atom_info%id_resname(j + 1) = str2id(s2s("__UNDEF__"))
175 38800 : atom_info%id_element(j + 1) = str2id(s2s("__UNDEF__"))
176 38800 : atom_info%resid(j + 1) = HUGE(0)
177 38800 : atom_info%atm_mass(j + 1) = HUGE(0.0_dp)
178 38800 : atom_info%atm_charge(j + 1) = -HUGE(0.0_dp)
179 38800 : atom_info%r(1, j + 1) = cp_unit_to_cp2k(atom_info%r(1, j + 1), "angstrom")
180 38800 : atom_info%r(2, j + 1) = cp_unit_to_cp2k(atom_info%r(2, j + 1), "angstrom")
181 38800 : atom_info%r(3, j + 1) = cp_unit_to_cp2k(atom_info%r(3, j + 1), "angstrom")
182 :
183 38826 : CALL parser_get_next_line(parser, 1, at_end=my_end)
184 : END DO
185 : ! Trigger error
186 26 : IF ((my_end) .AND. (j /= natom - MOD(natom, 2) + 1)) THEN
187 0 : IF (j /= natom) THEN
188 0 : CPABORT("Error while reading CRD file. Unexpected end of file.")
189 : END IF
190 26 : ELSE IF (MOD(natom, 2) /= 0) THEN
191 : ! In case let's handle the last atom
192 2 : j = natom
193 2 : 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 2 : atom_info%id_atmname(j) = str2id(s2s("__UNDEF__"))
197 2 : atom_info%id_molname(j) = str2id(s2s("__UNDEF__"))
198 2 : atom_info%id_resname(j) = str2id(s2s("__UNDEF__"))
199 2 : atom_info%id_element(j) = str2id(s2s("__UNDEF__"))
200 2 : atom_info%resid(j) = HUGE(0)
201 2 : atom_info%atm_mass(j) = HUGE(0.0_dp)
202 2 : atom_info%atm_charge(j) = -HUGE(0.0_dp)
203 2 : atom_info%r(1, j) = cp_unit_to_cp2k(atom_info%r(1, j), "angstrom")
204 2 : atom_info%r(2, j) = cp_unit_to_cp2k(atom_info%r(2, j), "angstrom")
205 2 : atom_info%r(3, j) = cp_unit_to_cp2k(atom_info%r(3, j), "angstrom")
206 :
207 2 : CALL parser_get_next_line(parser, 1, at_end=my_end)
208 : END IF
209 :
210 26 : IF (my_end) THEN
211 20 : CPWARN_IF(j /= natom, "No VELOCITY or BOX information found in CRD file.")
212 : ELSE
213 : ! Velocities
214 6 : CALL reallocate(velocity, 1, 3, 1, natom)
215 38604 : DO j = 1, natom - MOD(natom, 2), 2
216 38598 : IF (my_end) EXIT
217 38598 : READ (parser%input_line, *) velocity(1, j), velocity(2, j), velocity(3, j), &
218 77196 : velocity(1, j + 1), velocity(2, j + 1), velocity(3, j + 1)
219 :
220 38598 : velocity(1, j) = cp_unit_to_cp2k(velocity(1, j), "angstrom*ps^-1")
221 38598 : velocity(2, j) = cp_unit_to_cp2k(velocity(2, j), "angstrom*ps^-1")
222 38598 : velocity(3, j) = cp_unit_to_cp2k(velocity(3, j), "angstrom*ps^-1")
223 154392 : velocity(1:3, j) = velocity(1:3, j)*amber_conv_factor
224 :
225 38598 : velocity(1, j + 1) = cp_unit_to_cp2k(velocity(1, j + 1), "angstrom*ps^-1")
226 38598 : velocity(2, j + 1) = cp_unit_to_cp2k(velocity(2, j + 1), "angstrom*ps^-1")
227 38598 : velocity(3, j + 1) = cp_unit_to_cp2k(velocity(3, j + 1), "angstrom*ps^-1")
228 154392 : velocity(1:3, j + 1) = velocity(1:3, j + 1)*amber_conv_factor
229 :
230 38604 : CALL parser_get_next_line(parser, 1, at_end=my_end)
231 : END DO
232 6 : setup_velocities = .TRUE.
233 6 : IF ((my_end) .AND. (j /= natom - MOD(natom, 2) + 1)) THEN
234 0 : IF (j /= natom) THEN
235 : CALL cp_warn(__LOCATION__, &
236 : "No VELOCITY information found in CRD file. Ignoring BOX information. "// &
237 0 : "Please provide the BOX information directly from the main CP2K input! ")
238 : END IF
239 : setup_velocities = .FALSE.
240 6 : ELSE IF (MOD(natom, 2) /= 0) THEN
241 : ! In case let's handle the last atom
242 0 : j = natom
243 0 : READ (parser%input_line, *) velocity(1, j), velocity(2, j), velocity(3, j)
244 :
245 0 : velocity(1, j) = cp_unit_to_cp2k(velocity(1, j), "angstrom*ps^-1")
246 0 : velocity(2, j) = cp_unit_to_cp2k(velocity(2, j), "angstrom*ps^-1")
247 0 : velocity(3, j) = cp_unit_to_cp2k(velocity(3, j), "angstrom*ps^-1")
248 0 : velocity(1:3, j) = velocity(1:3, j)*amber_conv_factor
249 :
250 0 : CALL parser_get_next_line(parser, 1, at_end=my_end)
251 : END IF
252 : IF (setup_velocities) THEN
253 6 : velocity_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
254 : CALL section_velocity_val_set(velocity_section, velocity=velocity, &
255 6 : conv_factor=1.0_dp)
256 : END IF
257 6 : DEALLOCATE (velocity)
258 : END IF
259 26 : IF (my_end) THEN
260 20 : CPWARN_IF(j /= natom, "BOX information missing in CRD file.")
261 : ELSE
262 6 : IF (j /= natom) THEN
263 : CALL cp_warn(__LOCATION__, &
264 : "BOX information found in CRD file. They will be ignored. "// &
265 6 : "Please provide the BOX information directly from the main CP2K input!")
266 : END IF
267 : END IF
268 26 : CALL parser_release(parser)
269 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
270 26 : "PRINT%TOPOLOGY_INFO/CRD_INFO")
271 26 : CALL timestop(handle)
272 :
273 78 : 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 24 : 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 24 : NULLIFY (logger)
300 24 : CALL timeset(routineN, handle)
301 24 : logger => cp_get_default_logger()
302 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/AMBER_INFO", &
303 24 : extension=".subsysLog")
304 :
305 24 : atom_info => topology%atom_info
306 24 : 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 24 : atom_info=atom_info, conn_info=conn_info)
311 :
312 : ! Molnames have been internally generated
313 24 : topology%molname_generated = .TRUE.
314 :
315 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
316 24 : "PRINT%TOPOLOGY_INFO/AMBER_INFO")
317 24 : CALL timestop(handle)
318 24 : 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 40 : 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 80 : ALLOCATABLE, DIMENSION(:) :: isymbl, labres, raw_isymbl, 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 40 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iac, ib, ibh, icb, icbh, ico, icp, icph, &
391 40 : ict, icth, ip, iph, ipres, it, ith, &
392 40 : iwork, jb, jbh, jp, jph, jt, jth, kp, &
393 40 : kph, kt, kth, lp, lph
394 40 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: full_torsions
395 : LOGICAL :: check, valid_format
396 80 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: asol, bsol, cn1, cn2, phase, pk, pn, &
397 40 : req, rk, teq, tk
398 : TYPE(cp_parser_type) :: parser
399 :
400 40 : CALL timeset(routineN, handle)
401 40 : IF (output_unit > 0) WRITE (output_unit, '(/,A)') " AMBER_INFO| Reading Amber Topology File: "// &
402 3 : TRIM(filename)
403 40 : CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.TRUE.)
404 40 : valid_format = check_amber_8_std(parser, output_unit)
405 40 : IF (valid_format) THEN
406 1560 : DO WHILE (get_section_parmtop(parser, section, input_format))
407 40 : SELECT CASE (TRIM(section))
408 : CASE ("TITLE")
409 : ! Who cares about the title?
410 40 : CYCLE
411 : CASE ("POINTERS")
412 40 : 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 40 : natom = info(1)
416 40 : ntypes = info(2)
417 40 : nbonh = info(3)
418 40 : mbona = info(4)
419 40 : ntheth = info(5)
420 40 : mtheta = info(6)
421 40 : nphih = info(7)
422 40 : mphia = info(8)
423 40 : nhparm = info(9)
424 40 : nparm = info(10)
425 40 : nnb = info(11)
426 40 : nres = info(12)
427 40 : nbona = info(13)
428 40 : ntheta = info(14)
429 40 : nphia = info(15)
430 40 : numbnd = info(16)
431 40 : numang = info(17)
432 40 : nptra = info(18)
433 40 : natyp = info(19)
434 40 : nphb = info(20)
435 40 : ifpert = info(21)
436 40 : nbper = info(22)
437 40 : ngper = info(23)
438 40 : ndper = info(24)
439 40 : mbper = info(25)
440 40 : mgper = info(26)
441 40 : mdper = info(27)
442 40 : ifbox = info(28)
443 40 : nmxrs = info(29)
444 40 : ifcap = info(30)
445 40 : numextra = info(31)
446 :
447 : ! Print some info if requested
448 40 : IF (output_unit > 0) THEN
449 3 : WRITE (output_unit, '(A,/)') " AMBER_INFO| Information from AMBER topology file:"
450 : WRITE (output_unit, 1000) &
451 3 : natom, ntypes, nbonh, mbona, ntheth, mtheta, nphih, &
452 3 : mphia, nhparm, nparm, nnb, nres, nbona, ntheta, &
453 3 : nphia, numbnd, numang, nptra, natyp, nphb, ifbox, &
454 6 : nmxrs, ifcap, numextra
455 : END IF
456 :
457 : ! Allocate temporary arrays
458 40 : IF (do_connectivity) THEN
459 24 : check = PRESENT(atom_info) .AND. PRESENT(conn_info)
460 24 : CPASSERT(check)
461 24 : natom_prev = 0
462 24 : IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname)
463 : ! Allocate for extracting connectivity infos
464 72 : ALLOCATE (labres(nres))
465 72 : ALLOCATE (ipres(nres))
466 : END IF
467 40 : IF (do_forcefield) THEN
468 : ! Allocate for extracting forcefield infos
469 48 : ALLOCATE (iac(natom))
470 48 : ALLOCATE (ico(ntypes*ntypes))
471 48 : ALLOCATE (rk(numbnd))
472 32 : ALLOCATE (req(numbnd))
473 48 : ALLOCATE (tk(numang))
474 32 : ALLOCATE (teq(numang))
475 46 : ALLOCATE (pk(nptra))
476 30 : ALLOCATE (pn(nptra))
477 30 : ALLOCATE (phase(nptra))
478 48 : ALLOCATE (cn1(ntypes*(ntypes + 1)/2))
479 32 : ALLOCATE (cn2(ntypes*(ntypes + 1)/2))
480 32 : ALLOCATE (asol(ntypes*(ntypes + 1)/2))
481 32 : ALLOCATE (bsol(ntypes*(ntypes + 1)/2))
482 : END IF
483 : ! Always Allocate
484 120 : ALLOCATE (ibh(nbonh))
485 80 : ALLOCATE (jbh(nbonh))
486 80 : ALLOCATE (icbh(nbonh))
487 116 : ALLOCATE (ib(nbona))
488 76 : ALLOCATE (jb(nbona))
489 76 : ALLOCATE (icb(nbona))
490 120 : ALLOCATE (ith(ntheth))
491 80 : ALLOCATE (jth(ntheth))
492 80 : ALLOCATE (kth(ntheth))
493 80 : ALLOCATE (icth(ntheth))
494 116 : ALLOCATE (it(ntheta))
495 76 : ALLOCATE (jt(ntheta))
496 76 : ALLOCATE (kt(ntheta))
497 76 : ALLOCATE (ict(ntheta))
498 116 : ALLOCATE (iph(nphih))
499 76 : ALLOCATE (jph(nphih))
500 76 : ALLOCATE (kph(nphih))
501 76 : ALLOCATE (lph(nphih))
502 76 : ALLOCATE (icph(nphih))
503 116 : ALLOCATE (ip(nphia))
504 76 : ALLOCATE (jp(nphia))
505 76 : ALLOCATE (kp(nphia))
506 76 : ALLOCATE (lp(nphia))
507 76 : ALLOCATE (icp(nphia))
508 : CASE ("ATOM_NAME")
509 : ! Atom names are just ignored according the CP2K philosophy
510 40 : CYCLE
511 : CASE ("AMBER_ATOM_TYPE")
512 40 : IF (.NOT. do_connectivity) CYCLE
513 24 : CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom)
514 72 : ALLOCATE (strtmp_a(natom))
515 24 : CALL rd_amber_section(parser, section, strtmp_a, natom)
516 78896 : DO i = 1, natom
517 78896 : atom_info%id_atmname(natom_prev + i) = str2id(strtmp_a(i))
518 : END DO
519 24 : DEALLOCATE (strtmp_a)
520 : CASE ("CHARGE")
521 40 : IF (.NOT. do_connectivity) CYCLE
522 24 : CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom)
523 24 : CALL rd_amber_section(parser, section, atom_info%atm_charge(natom_prev + 1:), natom)
524 : ! Convert charges into atomic units
525 78896 : atom_info%atm_charge(natom_prev + 1:) = atom_info%atm_charge(natom_prev + 1:)/amber_conv_charge
526 : CASE ("MASS")
527 40 : IF (.NOT. do_connectivity) CYCLE
528 24 : CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom)
529 24 : CALL rd_amber_section(parser, section, atom_info%atm_mass(natom_prev + 1:), natom)
530 : CASE ("RESIDUE_LABEL")
531 40 : IF (.NOT. do_connectivity) CYCLE
532 24 : CALL reallocate(atom_info%id_resname, 1, natom_prev + natom)
533 24 : CALL rd_amber_section(parser, section, labres, nres)
534 : CASE ("RESIDUE_POINTER")
535 40 : IF (.NOT. do_connectivity) CYCLE
536 24 : CALL reallocate(atom_info%resid, 1, natom_prev + natom)
537 24 : CALL rd_amber_section(parser, section, ipres, nres)
538 : CASE ("ATOM_TYPE_INDEX")
539 40 : IF (.NOT. do_forcefield) CYCLE
540 16 : CALL rd_amber_section(parser, section, iac, natom)
541 : CASE ("NONBONDED_PARM_INDEX")
542 40 : IF (.NOT. do_forcefield) CYCLE
543 16 : CALL rd_amber_section(parser, section, ico, ntypes**2)
544 : CASE ("BOND_FORCE_CONSTANT")
545 40 : IF (.NOT. do_forcefield) CYCLE
546 16 : CALL rd_amber_section(parser, section, rk, numbnd)
547 : CASE ("BOND_EQUIL_VALUE")
548 40 : IF (.NOT. do_forcefield) CYCLE
549 16 : CALL rd_amber_section(parser, section, req, numbnd)
550 : CASE ("ANGLE_FORCE_CONSTANT")
551 40 : IF (.NOT. do_forcefield) CYCLE
552 16 : CALL rd_amber_section(parser, section, tk, numang)
553 : CASE ("ANGLE_EQUIL_VALUE")
554 40 : IF (.NOT. do_forcefield) CYCLE
555 16 : CALL rd_amber_section(parser, section, teq, numang)
556 : CASE ("DIHEDRAL_FORCE_CONSTANT")
557 40 : IF (.NOT. do_forcefield) CYCLE
558 16 : CALL rd_amber_section(parser, section, pk, nptra)
559 16 : IF (nptra <= 0) CYCLE
560 : ! Save raw values
561 14 : IF (ASSOCIATED(amb_info%raw_torsion_k)) DEALLOCATE (amb_info%raw_torsion_k)
562 380 : ALLOCATE (amb_info%raw_torsion_k(nptra), source=pk)
563 : CASE ("DIHEDRAL_PERIODICITY")
564 40 : IF (.NOT. do_forcefield) CYCLE
565 16 : CALL rd_amber_section(parser, section, pn, nptra)
566 16 : IF (nptra <= 0) CYCLE
567 : ! Save raw values
568 14 : IF (ASSOCIATED(amb_info%raw_torsion_m)) DEALLOCATE (amb_info%raw_torsion_m)
569 380 : ALLOCATE (amb_info%raw_torsion_m(nptra), source=pn)
570 : CASE ("DIHEDRAL_PHASE")
571 40 : IF (.NOT. do_forcefield) CYCLE
572 16 : CALL rd_amber_section(parser, section, phase, nptra)
573 16 : IF (nptra <= 0) CYCLE
574 : ! Save raw values
575 14 : IF (ASSOCIATED(amb_info%raw_torsion_phi0)) DEALLOCATE (amb_info%raw_torsion_phi0)
576 380 : ALLOCATE (amb_info%raw_torsion_phi0(nptra), source=phase)
577 : CASE ("LENNARD_JONES_ACOEF")
578 40 : IF (.NOT. do_forcefield) CYCLE
579 16 : CALL rd_amber_section(parser, section, cn1, ntypes*(ntypes + 1)/2)
580 : CASE ("LENNARD_JONES_BCOEF")
581 40 : IF (.NOT. do_forcefield) CYCLE
582 16 : CALL rd_amber_section(parser, section, cn2, ntypes*(ntypes + 1)/2)
583 : CASE ("HBOND_ACOEF")
584 40 : IF (.NOT. do_forcefield) CYCLE
585 16 : CALL rd_amber_section(parser, section, asol, nphb)
586 : CASE ("HBOND_BCOEF")
587 40 : IF (.NOT. do_forcefield) CYCLE
588 16 : 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 40 : CALL rd_amber_section(parser, section, ibh, jbh, icbh, nbonh)
592 : ! Conver to an atomic index
593 100064 : ibh(:) = ibh(:)/3 + 1
594 100064 : jbh(:) = jbh(:)/3 + 1
595 : CASE ("BONDS_WITHOUT_HYDROGEN")
596 : ! We always need to parse this information both for connectivity and forcefields
597 40 : CALL rd_amber_section(parser, section, ib, jb, icb, nbona)
598 : ! Conver to an atomic index
599 14058 : ib(:) = ib(:)/3 + 1
600 14058 : jb(:) = jb(:)/3 + 1
601 : CASE ("ANGLES_INC_HYDROGEN")
602 : ! We always need to parse this information both for connectivity and forcefields
603 40 : CALL rd_amber_section(parser, section, ith, jth, kth, icth, ntheth)
604 : ! Conver to an atomic index
605 72562 : ith(:) = ith(:)/3 + 1
606 72562 : jth(:) = jth(:)/3 + 1
607 72562 : kth(:) = kth(:)/3 + 1
608 : CASE ("ANGLES_WITHOUT_HYDROGEN")
609 : ! We always need to parse this information both for connectivity and forcefields
610 40 : CALL rd_amber_section(parser, section, it, jt, kt, ict, ntheta)
611 : ! Conver to an atomic index
612 18994 : it(:) = it(:)/3 + 1
613 18994 : jt(:) = jt(:)/3 + 1
614 18994 : kt(:) = kt(:)/3 + 1
615 : CASE ("DIHEDRALS_INC_HYDROGEN")
616 : ! We always need to parse this information both for connectivity and forcefields
617 40 : CALL rd_amber_section(parser, section, iph, jph, kph, lph, icph, nphih)
618 : ! Conver to an atomic index
619 56688 : iph(:) = iph(:)/3 + 1
620 56688 : jph(:) = jph(:)/3 + 1
621 56688 : kph(:) = ABS(kph(:))/3 + 1
622 56688 : lph(:) = ABS(lph(:))/3 + 1
623 : CASE ("DIHEDRALS_WITHOUT_HYDROGEN")
624 : ! We always need to parse this information both for connectivity and forcefields
625 40 : CALL rd_amber_section(parser, section, ip, jp, kp, lp, icp, nphia)
626 : ! Conver to an atomic index
627 45316 : ip(:) = ip(:)/3 + 1
628 45316 : jp(:) = jp(:)/3 + 1
629 45316 : kp(:) = ABS(kp(:))/3 + 1
630 46852 : 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 40 : IF (do_forcefield .AND. (nphih + nphia > 0)) THEN
637 14 : IF (ASSOCIATED(amb_info%raw_torsion_id)) DEALLOCATE (amb_info%raw_torsion_id)
638 42 : ALLOCATE (amb_info%raw_torsion_id(5, nphih + nphia))
639 28132 : DO i = 1, nphih
640 28118 : amb_info%raw_torsion_id(1, i) = iph(i)
641 28118 : amb_info%raw_torsion_id(2, i) = jph(i)
642 28118 : amb_info%raw_torsion_id(3, i) = kph(i)
643 28118 : amb_info%raw_torsion_id(4, i) = lph(i)
644 28132 : amb_info%raw_torsion_id(5, i) = icph(i)
645 : END DO
646 22356 : DO i = 1, nphia
647 22342 : amb_info%raw_torsion_id(1, nphih + i) = ip(i)
648 22342 : amb_info%raw_torsion_id(2, nphih + i) = jp(i)
649 22342 : amb_info%raw_torsion_id(3, nphih + i) = kp(i)
650 22342 : amb_info%raw_torsion_id(4, nphih + i) = lp(i)
651 22356 : 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 40 : IF (do_connectivity) THEN
658 24 : CALL timeset(TRIM(routineN)//"_connectivity", handle2)
659 : ! ----------------------------------------------------------
660 : ! Conform Amber Names with CHARMM convention (kind<->charge)
661 : ! ----------------------------------------------------------
662 72 : ALLOCATE (isymbl(natom))
663 72 : ALLOCATE (iwork(natom))
664 :
665 78896 : DO i = 1, SIZE(isymbl)
666 78896 : isymbl(i) = id2str(atom_info%id_atmname(natom_prev + i))
667 : END DO
668 :
669 : ! Sort atom names + charges and identify unique types
670 24 : CALL sort(isymbl, natom, iwork)
671 48 : ALLOCATE (raw_isymbl(natom))
672 78896 : raw_isymbl(:) = isymbl
673 :
674 24 : istart = 1
675 78872 : DO i = 2, natom
676 78872 : IF (TRIM(isymbl(i)) /= TRIM(isymbl(istart))) THEN
677 244 : CALL conform_atom_type_low(isymbl, iwork, i, istart, atom_info%atm_charge(natom_prev + 1:))
678 244 : istart = i
679 : END IF
680 : END DO
681 24 : CALL conform_atom_type_low(isymbl, iwork, i, istart, atom_info%atm_charge(natom_prev + 1:))
682 :
683 24 : CALL resolve_atom_type_collisions(isymbl, raw_isymbl)
684 :
685 : ! Copy back the modified and conformed atom types
686 78896 : DO i = 1, natom
687 78896 : atom_info%id_atmname(natom_prev + iwork(i)) = str2id(s2s(isymbl(i)))
688 : END DO
689 :
690 : ! -----------------------------------------------------------
691 : ! Fill residue_name and residue_id information before exiting
692 : ! -----------------------------------------------------------
693 22734 : DO i = 1, nres - 1
694 123798 : atom_info%id_resname(natom_prev + ipres(i):natom_prev + ipres(i + 1)) = str2id(s2s(labres(i)))
695 123822 : atom_info%resid(natom_prev + ipres(i):natom_prev + ipres(i + 1)) = i
696 : END DO
697 518 : atom_info%id_resname(natom_prev + ipres(i):natom_prev + natom) = str2id(s2s(labres(i)))
698 518 : atom_info%resid(natom_prev + ipres(i):natom_prev + natom) = i
699 :
700 : ! Deallocate when extracting connectivity infos
701 24 : DEALLOCATE (iwork)
702 24 : DEALLOCATE (isymbl)
703 24 : DEALLOCATE (raw_isymbl)
704 24 : DEALLOCATE (labres)
705 24 : DEALLOCATE (ipres)
706 :
707 : ! ----------------------------------------------------------
708 : ! Copy connectivity
709 : ! ----------------------------------------------------------
710 : ! BONDS
711 24 : nbond_prev = 0
712 24 : IF (ASSOCIATED(conn_info%bond_a)) nbond_prev = SIZE(conn_info%bond_a)
713 :
714 24 : CALL reallocate(conn_info%bond_a, 1, nbond_prev + nbonh + nbona)
715 24 : CALL reallocate(conn_info%bond_b, 1, nbond_prev + nbonh + nbona)
716 50096 : DO i = 1, nbonh
717 50072 : index_now = nbond_prev + i
718 50072 : conn_info%bond_a(index_now) = natom_prev + ibh(i)
719 50096 : conn_info%bond_b(index_now) = natom_prev + jbh(i)
720 : END DO
721 7162 : DO i = 1, nbona
722 7138 : index_now = nbond_prev + i + nbonh
723 7138 : conn_info%bond_a(index_now) = natom_prev + ib(i)
724 7162 : conn_info%bond_b(index_now) = natom_prev + jb(i)
725 : END DO
726 :
727 : ! ANGLES
728 24 : ntheta_prev = 0
729 24 : IF (ASSOCIATED(conn_info%theta_a)) ntheta_prev = SIZE(conn_info%theta_a)
730 :
731 24 : CALL reallocate(conn_info%theta_a, 1, ntheta_prev + ntheth + ntheta)
732 24 : CALL reallocate(conn_info%theta_b, 1, ntheta_prev + ntheth + ntheta)
733 24 : CALL reallocate(conn_info%theta_c, 1, ntheta_prev + ntheth + ntheta)
734 36406 : DO i = 1, ntheth
735 36382 : index_now = ntheta_prev + i
736 36382 : conn_info%theta_a(index_now) = natom_prev + ith(i)
737 36382 : conn_info%theta_b(index_now) = natom_prev + jth(i)
738 36406 : conn_info%theta_c(index_now) = natom_prev + kth(i)
739 : END DO
740 9692 : DO i = 1, ntheta
741 9668 : index_now = ntheta_prev + i + ntheth
742 9668 : conn_info%theta_a(index_now) = natom_prev + it(i)
743 9668 : conn_info%theta_b(index_now) = natom_prev + jt(i)
744 9692 : conn_info%theta_c(index_now) = natom_prev + kt(i)
745 : END DO
746 :
747 : ! TORSIONS
748 : ! For torsions we need to find out the unique torsions
749 : ! defined in the amber parmtop
750 24 : nphi_prev = 0
751 24 : IF (ASSOCIATED(conn_info%phi_a)) nphi_prev = SIZE(conn_info%phi_a)
752 :
753 24 : CALL reallocate(conn_info%phi_a, 1, nphi_prev + nphih + nphia)
754 24 : CALL reallocate(conn_info%phi_b, 1, nphi_prev + nphih + nphia)
755 24 : CALL reallocate(conn_info%phi_c, 1, nphi_prev + nphih + nphia)
756 24 : CALL reallocate(conn_info%phi_d, 1, nphi_prev + nphih + nphia)
757 :
758 24 : IF (nphih + nphia /= 0) THEN
759 66 : ALLOCATE (full_torsions(4, nphih + nphia))
760 66 : ALLOCATE (iwork(nphih + nphia))
761 :
762 28552 : DO i = 1, nphih
763 28530 : full_torsions(1, i) = iph(i)
764 28530 : full_torsions(2, i) = jph(i)
765 28530 : full_torsions(3, i) = kph(i)
766 28552 : full_torsions(4, i) = lph(i)
767 : END DO
768 22956 : DO i = 1, nphia
769 22934 : full_torsions(1, nphih + i) = ip(i)
770 22934 : full_torsions(2, nphih + i) = jp(i)
771 22934 : full_torsions(3, nphih + i) = kp(i)
772 22956 : full_torsions(4, nphih + i) = lp(i)
773 : END DO
774 22 : CALL sort(full_torsions, 1, nphih + nphia, 1, 4, iwork)
775 :
776 22 : unique_torsions = nphi_prev + 1
777 22 : conn_info%phi_a(unique_torsions) = natom_prev + full_torsions(1, 1)
778 22 : conn_info%phi_b(unique_torsions) = natom_prev + full_torsions(2, 1)
779 22 : conn_info%phi_c(unique_torsions) = natom_prev + full_torsions(3, 1)
780 22 : conn_info%phi_d(unique_torsions) = natom_prev + full_torsions(4, 1)
781 51464 : DO i = 2, nphih + nphia
782 : IF ((full_torsions(1, i) /= full_torsions(1, i - 1)) .OR. &
783 : (full_torsions(2, i) /= full_torsions(2, i - 1)) .OR. &
784 51442 : (full_torsions(3, i) /= full_torsions(3, i - 1)) .OR. &
785 22 : (full_torsions(4, i) /= full_torsions(4, i - 1))) THEN
786 37652 : unique_torsions = unique_torsions + 1
787 37652 : conn_info%phi_a(unique_torsions) = natom_prev + full_torsions(1, i)
788 37652 : conn_info%phi_b(unique_torsions) = natom_prev + full_torsions(2, i)
789 37652 : conn_info%phi_c(unique_torsions) = natom_prev + full_torsions(3, i)
790 37652 : conn_info%phi_d(unique_torsions) = natom_prev + full_torsions(4, i)
791 : END IF
792 : END DO
793 22 : CALL reallocate(conn_info%phi_a, 1, unique_torsions)
794 22 : CALL reallocate(conn_info%phi_b, 1, unique_torsions)
795 22 : CALL reallocate(conn_info%phi_c, 1, unique_torsions)
796 22 : CALL reallocate(conn_info%phi_d, 1, unique_torsions)
797 :
798 22 : DEALLOCATE (full_torsions)
799 22 : DEALLOCATE (iwork)
800 : END IF
801 : ! IMPROPERS
802 24 : CALL reallocate(conn_info%impr_a, 1, 0)
803 24 : CALL reallocate(conn_info%impr_b, 1, 0)
804 24 : CALL reallocate(conn_info%impr_c, 1, 0)
805 24 : CALL reallocate(conn_info%impr_d, 1, 0)
806 :
807 : ! ----------------------------------------------------------
808 : ! Generate molecule names
809 : ! ----------------------------------------------------------
810 24 : CALL reallocate(atom_info%id_molname, 1, natom_prev + natom)
811 78896 : atom_info%id_molname(natom_prev + 1:natom_prev + natom) = str2id(s2s("__UNDEF__"))
812 : CALL topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, &
813 24 : atom_info%id_molname(natom_prev + 1:natom_prev + natom))
814 48 : CALL timestop(handle2)
815 : END IF
816 :
817 : ! Extracts force fields info from the AMBER topology file
818 40 : IF (do_forcefield) THEN
819 16 : CALL timeset(TRIM(routineN)//"_forcefield", handle2)
820 : ! ----------------------------------------------------------
821 : ! Force Fields informations related to bonds
822 : ! ----------------------------------------------------------
823 16 : CALL reallocate(amb_info%bond_a, 1, buffer_size)
824 16 : CALL reallocate(amb_info%bond_b, 1, buffer_size)
825 16 : CALL reallocate(amb_info%bond_k, 1, buffer_size)
826 16 : CALL reallocate(amb_info%bond_r0, 1, buffer_size)
827 16 : nsize = 0
828 : ! Bonds containing hydrogens
829 : CALL post_process_bonds_info(amb_info%bond_a, amb_info%bond_b, &
830 : amb_info%bond_k, amb_info%bond_r0, particle_set, nsize, &
831 16 : nbonh, ibh, jbh, icbh, rk, req)
832 : ! Bonds non-containing hydrogens
833 : CALL post_process_bonds_info(amb_info%bond_a, amb_info%bond_b, &
834 : amb_info%bond_k, amb_info%bond_r0, particle_set, nsize, &
835 16 : nbona, ib, jb, icb, rk, req)
836 : ! Shrink arrays size to the minimal request
837 16 : CALL reallocate(amb_info%bond_a, 1, nsize)
838 16 : CALL reallocate(amb_info%bond_b, 1, nsize)
839 16 : CALL reallocate(amb_info%bond_k, 1, nsize)
840 16 : CALL reallocate(amb_info%bond_r0, 1, nsize)
841 :
842 : ! ----------------------------------------------------------
843 : ! Force Fields informations related to bends
844 : ! ----------------------------------------------------------
845 16 : CALL reallocate(amb_info%bend_a, 1, buffer_size)
846 16 : CALL reallocate(amb_info%bend_b, 1, buffer_size)
847 16 : CALL reallocate(amb_info%bend_c, 1, buffer_size)
848 16 : CALL reallocate(amb_info%bend_k, 1, buffer_size)
849 16 : CALL reallocate(amb_info%bend_theta0, 1, buffer_size)
850 16 : nsize = 0
851 : ! Bends containing hydrogens
852 : CALL post_process_bends_info(amb_info%bend_a, amb_info%bend_b, &
853 : amb_info%bend_c, amb_info%bend_k, amb_info%bend_theta0, &
854 16 : particle_set, nsize, ntheth, ith, jth, kth, icth, tk, teq)
855 : ! Bends non-containing hydrogens
856 : CALL post_process_bends_info(amb_info%bend_a, amb_info%bend_b, &
857 : amb_info%bend_c, amb_info%bend_k, amb_info%bend_theta0, &
858 16 : particle_set, nsize, ntheta, it, jt, kt, ict, tk, teq)
859 : ! Shrink arrays size to the minimal request
860 16 : CALL reallocate(amb_info%bend_a, 1, nsize)
861 16 : CALL reallocate(amb_info%bend_b, 1, nsize)
862 16 : CALL reallocate(amb_info%bend_c, 1, nsize)
863 16 : CALL reallocate(amb_info%bend_k, 1, nsize)
864 16 : CALL reallocate(amb_info%bend_theta0, 1, nsize)
865 :
866 : ! ----------------------------------------------------------
867 : ! Force Fields informations related to torsions
868 : ! in amb_info%phi0 we store PHI0
869 : ! ----------------------------------------------------------
870 :
871 16 : CALL reallocate(amb_info%torsion_a, 1, buffer_size)
872 16 : CALL reallocate(amb_info%torsion_b, 1, buffer_size)
873 16 : CALL reallocate(amb_info%torsion_c, 1, buffer_size)
874 16 : CALL reallocate(amb_info%torsion_d, 1, buffer_size)
875 16 : CALL reallocate(amb_info%torsion_k, 1, buffer_size)
876 16 : CALL reallocate(amb_info%torsion_m, 1, buffer_size)
877 16 : CALL reallocate(amb_info%torsion_phi0, 1, buffer_size)
878 16 : nsize = 0
879 : ! Torsions 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 16 : nphih, iph, jph, kph, lph, icph, pk, pn, phase)
884 : ! Torsions non-containing hydrogens
885 : CALL post_process_torsions_info(amb_info%torsion_a, amb_info%torsion_b, &
886 : amb_info%torsion_c, amb_info%torsion_d, amb_info%torsion_k, &
887 : amb_info%torsion_m, amb_info%torsion_phi0, particle_set, nsize, &
888 16 : nphia, ip, jp, kp, lp, icp, pk, pn, phase)
889 : ! Shrink arrays size to the minimal request
890 16 : CALL reallocate(amb_info%torsion_a, 1, nsize)
891 16 : CALL reallocate(amb_info%torsion_b, 1, nsize)
892 16 : CALL reallocate(amb_info%torsion_c, 1, nsize)
893 16 : CALL reallocate(amb_info%torsion_d, 1, nsize)
894 16 : CALL reallocate(amb_info%torsion_k, 1, nsize)
895 16 : CALL reallocate(amb_info%torsion_m, 1, nsize)
896 16 : CALL reallocate(amb_info%torsion_phi0, 1, nsize)
897 :
898 : ! Sort dihedral metadata for faster lookup
899 16 : IF (nphih + nphia /= 0) THEN
900 42 : ALLOCATE (iwork(nphih + nphia))
901 14 : CALL sort(amb_info%raw_torsion_id, 1, nphih + nphia, 1, 5, iwork)
902 14 : DEALLOCATE (iwork)
903 : END IF
904 :
905 : ! ----------------------------------------------------------
906 : ! Post process of LJ parameters
907 : ! ----------------------------------------------------------
908 16 : CALL reallocate(amb_info%nonbond_a, 1, buffer_size)
909 16 : CALL reallocate(amb_info%nonbond_eps, 1, buffer_size)
910 16 : CALL reallocate(amb_info%nonbond_rmin2, 1, buffer_size)
911 :
912 16 : nsize = 0
913 : CALL post_process_LJ_info(amb_info%nonbond_a, amb_info%nonbond_eps, &
914 : amb_info%nonbond_rmin2, particle_set, ntypes, nsize, iac, ico, &
915 16 : cn1, cn2, natom)
916 :
917 : ! Shrink arrays size to the minimal request
918 16 : CALL reallocate(amb_info%nonbond_a, 1, nsize)
919 16 : CALL reallocate(amb_info%nonbond_eps, 1, nsize)
920 16 : CALL reallocate(amb_info%nonbond_rmin2, 1, nsize)
921 :
922 : ! Deallocate at the end of the dirty job
923 16 : DEALLOCATE (iac)
924 16 : DEALLOCATE (ico)
925 16 : DEALLOCATE (rk)
926 16 : DEALLOCATE (req)
927 16 : DEALLOCATE (tk)
928 16 : DEALLOCATE (teq)
929 16 : DEALLOCATE (pk)
930 16 : DEALLOCATE (pn)
931 16 : DEALLOCATE (phase)
932 16 : DEALLOCATE (cn1)
933 16 : DEALLOCATE (cn2)
934 16 : DEALLOCATE (asol)
935 16 : DEALLOCATE (bsol)
936 16 : CALL timestop(handle2)
937 : END IF
938 : ! Always Deallocate
939 40 : DEALLOCATE (ibh)
940 40 : DEALLOCATE (jbh)
941 40 : DEALLOCATE (icbh)
942 40 : DEALLOCATE (ib)
943 40 : DEALLOCATE (jb)
944 40 : DEALLOCATE (icb)
945 40 : DEALLOCATE (ith)
946 40 : DEALLOCATE (jth)
947 40 : DEALLOCATE (kth)
948 40 : DEALLOCATE (icth)
949 40 : DEALLOCATE (it)
950 40 : DEALLOCATE (jt)
951 40 : DEALLOCATE (kt)
952 40 : DEALLOCATE (ict)
953 40 : DEALLOCATE (iph)
954 40 : DEALLOCATE (jph)
955 40 : DEALLOCATE (kph)
956 40 : DEALLOCATE (lph)
957 40 : DEALLOCATE (icph)
958 40 : DEALLOCATE (ip)
959 40 : DEALLOCATE (jp)
960 40 : DEALLOCATE (kp)
961 40 : DEALLOCATE (lp)
962 40 : DEALLOCATE (icp)
963 40 : CALL parser_release(parser)
964 40 : CALL timestop(handle)
965 40 : RETURN
966 : ! Output info Format
967 : 1000 FORMAT(T2, &
968 : /' NATOM = ', i7, ' NTYPES = ', i7, ' NBONH = ', i7, ' MBONA = ', i7, &
969 : /' NTHETH = ', i7, ' MTHETA = ', i7, ' NPHIH = ', i7, ' MPHIA = ', i7, &
970 : /' NHPARM = ', i7, ' NPARM = ', i7, ' NNB = ', i7, ' NRES = ', i7, &
971 : /' NBONA = ', i7, ' NTHETA = ', i7, ' NPHIA = ', i7, ' NUMBND = ', i7, &
972 : /' NUMANG = ', i7, ' NPTRA = ', i7, ' NATYP = ', i7, ' NPHB = ', i7, &
973 : /' IFBOX = ', i7, ' NMXRS = ', i7, ' IFCAP = ', i7, ' NEXTRA = ', i7,/)
974 160 : END SUBROUTINE rdparm_amber_8
975 :
976 : ! **************************************************************************************************
977 : !> \brief Low level routine to identify and rename unique atom types
978 : !> \param isymbl ...
979 : !> \param iwork ...
980 : !> \param i ...
981 : !> \param istart ...
982 : !> \param charges ...
983 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
984 : ! **************************************************************************************************
985 268 : SUBROUTINE conform_atom_type_low(isymbl, iwork, i, istart, charges)
986 : CHARACTER(LEN=default_string_length), DIMENSION(:) :: isymbl
987 : INTEGER, DIMENSION(:) :: iwork
988 : INTEGER, INTENT(IN) :: i
989 : INTEGER, INTENT(INOUT) :: istart
990 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges
991 :
992 : CHARACTER(len=*), PARAMETER :: routineN = 'conform_atom_type_low'
993 :
994 : INTEGER :: counter, gind, handle, iend, ind, isize, &
995 : j, k, kend, kstart
996 268 : INTEGER, DIMENSION(:), POINTER :: cindx, lindx
997 : REAL(KIND=dp) :: ctmp
998 268 : REAL(KIND=dp), DIMENSION(:), POINTER :: cwork
999 :
1000 268 : CALL timeset(routineN, handle)
1001 268 : iend = i - 1
1002 268 : isize = iend - istart + 1
1003 804 : ALLOCATE (cwork(isize))
1004 804 : ALLOCATE (lindx(isize))
1005 536 : ALLOCATE (cindx(isize))
1006 268 : ind = 0
1007 79140 : DO k = istart, iend
1008 78872 : ind = ind + 1
1009 78872 : cwork(ind) = charges(iwork(k))
1010 79140 : lindx(ind) = k
1011 : END DO
1012 268 : CALL sort(cwork, isize, cindx)
1013 :
1014 268 : ctmp = cwork(1)
1015 268 : counter = 1
1016 78872 : DO k = 2, isize
1017 78872 : IF (cwork(k) /= ctmp) THEN
1018 1414 : counter = counter + 1
1019 1414 : ctmp = cwork(k)
1020 : END IF
1021 : END DO
1022 268 : IF (counter /= 1) THEN
1023 154 : counter = 1
1024 154 : kstart = 1
1025 154 : ctmp = cwork(1)
1026 12778 : DO k = 2, isize
1027 12778 : IF (cwork(k) /= ctmp) THEN
1028 : kend = k - 1
1029 12364 : DO j = kstart, kend
1030 10950 : gind = lindx(cindx(j))
1031 10950 : IF (LEN_TRIM(isymbl(gind)) + LEN_TRIM(ADJUSTL(cp_to_string(counter))) > &
1032 10950 : LEN(isymbl(gind))) CPABORT("AMBER charge kind name exceeds string length")
1033 12364 : isymbl(gind) = TRIM(isymbl(gind))//ADJUSTL(cp_to_string(counter))
1034 : END DO
1035 1414 : counter = counter + 1
1036 1414 : ctmp = cwork(k)
1037 1414 : kstart = k
1038 : END IF
1039 : END DO
1040 : kend = k - 1
1041 1982 : DO j = kstart, kend
1042 1828 : gind = lindx(cindx(j))
1043 1828 : IF (LEN_TRIM(isymbl(gind)) + LEN_TRIM(ADJUSTL(cp_to_string(counter))) > &
1044 1828 : LEN(isymbl(gind))) CPABORT("AMBER charge kind name exceeds string length")
1045 1982 : isymbl(gind) = TRIM(isymbl(gind))//ADJUSTL(cp_to_string(counter))
1046 : END DO
1047 : END IF
1048 268 : DEALLOCATE (cwork)
1049 268 : DEALLOCATE (lindx)
1050 268 : DEALLOCATE (cindx)
1051 268 : CALL timestop(handle)
1052 268 : END SUBROUTINE conform_atom_type_low
1053 :
1054 : ! **************************************************************************************************
1055 : !> \brief Resolve ambiguities introduced by the legacy AMBER charge-split
1056 : !> atom-kind naming scheme while preserving non-colliding names.
1057 : !> \param isymbl Conformed atom-kind names
1058 : !> \param raw_isymbl Original AMBER atom types
1059 : ! **************************************************************************************************
1060 24 : SUBROUTINE resolve_atom_type_collisions(isymbl, raw_isymbl)
1061 : CHARACTER(LEN=default_string_length), &
1062 : DIMENSION(:), INTENT(INOUT) :: isymbl
1063 : CHARACTER(LEN=default_string_length), &
1064 : DIMENSION(:), INTENT(IN) :: raw_isymbl
1065 :
1066 : CHARACTER(len=*), PARAMETER :: routineN = 'resolve_atom_type_collisions'
1067 :
1068 : CHARACTER(LEN=default_string_length) :: base_name, candidate, suffix
1069 : CHARACTER(LEN=default_string_length), &
1070 24 : ALLOCATABLE, DIMENSION(:) :: sorted_labels
1071 : INTEGER :: finish, handle, idx, jdx, keeper, n, &
1072 : nsep, p, q, raw_len, start
1073 24 : INTEGER, ALLOCATABLE, DIMENSION(:) :: order
1074 : LOGICAL :: collision, seen_type
1075 :
1076 24 : CALL timeset(routineN, handle)
1077 :
1078 24 : n = SIZE(isymbl)
1079 24 : CPASSERT(SIZE(raw_isymbl) == n)
1080 :
1081 72 : ALLOCATE (sorted_labels(n))
1082 72 : ALLOCATE (order(n))
1083 :
1084 : ! Sort a copy of the legacy-generated labels. Equal labels then
1085 : ! form contiguous groups, while order maps back to isymbl/raw_isymbl.
1086 78896 : sorted_labels(:) = isymbl
1087 24 : CALL sort(sorted_labels, n, order)
1088 :
1089 24 : start = 1
1090 1680 : DO WHILE (start <= n)
1091 : finish = start
1092 :
1093 78872 : DO WHILE (finish < n)
1094 78848 : IF (sorted_labels(finish + 1) /= sorted_labels(start)) EXIT
1095 1656 : finish = finish + 1
1096 : END DO
1097 :
1098 1656 : IF (finish > start) THEN
1099 1372 : keeper = order(start)
1100 1372 : collision = .FALSE.
1101 :
1102 : ! Repeated atoms of the same original AMBER type are normal.
1103 : ! A collision exists only when one generated label represents
1104 : ! more than one original AMBER atom type.
1105 78392 : DO p = start + 1, finish
1106 77046 : idx = order(p)
1107 78392 : IF (raw_isymbl(idx) /= raw_isymbl(keeper)) THEN
1108 : collision = .TRUE.
1109 : EXIT
1110 : END IF
1111 : END DO
1112 :
1113 1372 : IF (collision) THEN
1114 :
1115 : ! Preserve an unsplit original AMBER atom type when its name
1116 : ! exactly matches the colliding legacy label. If no such type
1117 : ! exists, disambiguate all identities in the collision group.
1118 734 : keeper = 0
1119 734 : DO p = start, finish
1120 716 : idx = order(p)
1121 734 : IF (raw_isymbl(idx) == sorted_labels(start)) THEN
1122 : keeper = idx
1123 : EXIT
1124 : END IF
1125 : END DO
1126 :
1127 744 : DO p = start, finish
1128 718 : idx = order(p)
1129 :
1130 718 : IF (keeper /= 0) THEN
1131 64 : IF (raw_isymbl(idx) == raw_isymbl(keeper)) CYCLE
1132 : END IF
1133 :
1134 : ! Handle each conflicting original AMBER type only once.
1135 3660 : seen_type = .FALSE.
1136 3660 : DO q = start, p - 1
1137 3616 : jdx = order(q)
1138 3660 : IF (raw_isymbl(jdx) == raw_isymbl(idx)) THEN
1139 : seen_type = .TRUE.
1140 : EXIT
1141 : END IF
1142 : END DO
1143 708 : IF (seen_type) CYCLE
1144 :
1145 44 : base_name = TRIM(raw_isymbl(idx))
1146 44 : raw_len = LEN_TRIM(base_name)
1147 :
1148 : ! A non-keeper colliding identity must have acquired a
1149 : ! numeric charge-class suffix through conform_atom_type_low.
1150 44 : IF (LEN_TRIM(sorted_labels(start)) <= raw_len) THEN
1151 0 : CPABORT("Unable to disambiguate AMBER atom kind name")
1152 : END IF
1153 :
1154 44 : suffix = ADJUSTL(sorted_labels(start) (raw_len + 1:))
1155 44 : nsep = 1
1156 :
1157 : ! Prefer TYPE_1. If that name already exists, use
1158 : ! TYPE__1, TYPE___1, ... until an unused name is found.
1159 0 : DO
1160 44 : IF (raw_len + nsep + LEN_TRIM(suffix) > LEN(candidate)) THEN
1161 0 : CPABORT("AMBER charge kind name exceeds string length")
1162 : END IF
1163 :
1164 88 : candidate = TRIM(base_name)//REPEAT("_", nsep)//TRIM(suffix)
1165 :
1166 540450 : IF (.NOT. ANY(isymbl == candidate)) EXIT
1167 0 : nsep = nsep + 1
1168 : END DO
1169 :
1170 : ! Rename all atoms belonging to this conflicting identity.
1171 1442 : DO q = start, finish
1172 1372 : jdx = order(q)
1173 2090 : IF (raw_isymbl(jdx) == raw_isymbl(idx)) THEN
1174 708 : isymbl(jdx) = candidate
1175 : END IF
1176 : END DO
1177 :
1178 : END DO
1179 : END IF
1180 : END IF
1181 :
1182 1656 : start = finish + 1
1183 : END DO
1184 :
1185 24 : DEALLOCATE (sorted_labels)
1186 24 : DEALLOCATE (order)
1187 :
1188 24 : CALL timestop(handle)
1189 24 : END SUBROUTINE resolve_atom_type_collisions
1190 :
1191 : ! **************************************************************************************************
1192 : !> \brief Set of Low level subroutines reading section for parmtop
1193 : !> reading 1 array of integers of length dim
1194 : !> \param parser ...
1195 : !> \param section ...
1196 : !> \param array1 ...
1197 : !> \param dim ...
1198 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1199 : ! **************************************************************************************************
1200 96 : SUBROUTINE rd_amber_section_i1(parser, section, array1, dim)
1201 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1202 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1203 : INTEGER, DIMENSION(:) :: array1
1204 : INTEGER, INTENT(IN) :: dim
1205 :
1206 : INTEGER :: i
1207 : LOGICAL :: my_end
1208 :
1209 96 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1210 96 : i = 1
1211 104626 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1212 104530 : IF (parser_test_next_token(parser) == "EOL") THEN
1213 114948 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1214 : END IF
1215 104530 : IF (my_end) EXIT
1216 104530 : CALL parser_get_object(parser, array1(i))
1217 104530 : i = i + 1
1218 : END DO
1219 : ! Trigger end of file aborting
1220 96 : IF (my_end .AND. (i <= dim)) THEN
1221 : CALL cp_abort(__LOCATION__, &
1222 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1223 : END IF
1224 96 : END SUBROUTINE rd_amber_section_i1
1225 :
1226 : ! **************************************************************************************************
1227 : !> \brief Set of Low level subroutines reading section for parmtop
1228 : !> reading 3 arrays of integers of length dim
1229 : !> \param parser ...
1230 : !> \param section ...
1231 : !> \param array1 ...
1232 : !> \param array2 ...
1233 : !> \param array3 ...
1234 : !> \param dim ...
1235 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1236 : ! **************************************************************************************************
1237 80 : SUBROUTINE rd_amber_section_i3(parser, section, array1, array2, array3, dim)
1238 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1239 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1240 : INTEGER, DIMENSION(:) :: array1, array2, array3
1241 : INTEGER, INTENT(IN) :: dim
1242 :
1243 : INTEGER :: i
1244 : LOGICAL :: my_end
1245 :
1246 80 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1247 80 : i = 1
1248 114122 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1249 : !array1
1250 114042 : IF (parser_test_next_token(parser) == "EOL") THEN
1251 125400 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1252 : END IF
1253 114042 : IF (my_end) EXIT
1254 114042 : CALL parser_get_object(parser, array1(i))
1255 : !array2
1256 114042 : IF (parser_test_next_token(parser) == "EOL") THEN
1257 125470 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1258 : END IF
1259 114042 : IF (my_end) EXIT
1260 114042 : CALL parser_get_object(parser, array2(i))
1261 : !array3
1262 114042 : IF (parser_test_next_token(parser) == "EOL") THEN
1263 125428 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1264 : END IF
1265 114042 : IF (my_end) EXIT
1266 114042 : CALL parser_get_object(parser, array3(i))
1267 114042 : i = i + 1
1268 : END DO
1269 : ! Trigger end of file aborting
1270 80 : IF (my_end .AND. (i <= dim)) THEN
1271 : CALL cp_abort(__LOCATION__, &
1272 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1273 : END IF
1274 80 : END SUBROUTINE rd_amber_section_i3
1275 :
1276 : ! **************************************************************************************************
1277 : !> \brief Set of Low level subroutines reading section for parmtop
1278 : !> reading 4 arrays of integers of length dim
1279 : !> \param parser ...
1280 : !> \param section ...
1281 : !> \param array1 ...
1282 : !> \param array2 ...
1283 : !> \param array3 ...
1284 : !> \param array4 ...
1285 : !> \param dim ...
1286 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1287 : ! **************************************************************************************************
1288 80 : SUBROUTINE rd_amber_section_i4(parser, section, array1, array2, array3, array4, dim)
1289 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1290 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1291 : INTEGER, DIMENSION(:) :: array1, array2, array3, array4
1292 : INTEGER, INTENT(IN) :: dim
1293 :
1294 : INTEGER :: i
1295 : LOGICAL :: my_end
1296 :
1297 80 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1298 80 : i = 1
1299 91556 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1300 : !array1
1301 91476 : IF (parser_test_next_token(parser) == "EOL") THEN
1302 109730 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1303 : END IF
1304 91476 : IF (my_end) EXIT
1305 91476 : CALL parser_get_object(parser, array1(i))
1306 : !array2
1307 91476 : IF (parser_test_next_token(parser) == "EOL") THEN
1308 91476 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1309 : END IF
1310 91476 : IF (my_end) EXIT
1311 91476 : CALL parser_get_object(parser, array2(i))
1312 : !array3
1313 91476 : IF (parser_test_next_token(parser) == "EOL") THEN
1314 109766 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1315 : END IF
1316 91476 : IF (my_end) EXIT
1317 91476 : CALL parser_get_object(parser, array3(i))
1318 : !array4
1319 91476 : IF (parser_test_next_token(parser) == "EOL") THEN
1320 91476 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1321 : END IF
1322 91476 : IF (my_end) EXIT
1323 91476 : CALL parser_get_object(parser, array4(i))
1324 91476 : i = i + 1
1325 : END DO
1326 : ! Trigger end of file aborting
1327 80 : IF (my_end .AND. (i <= dim)) THEN
1328 : CALL cp_abort(__LOCATION__, &
1329 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1330 : END IF
1331 80 : END SUBROUTINE rd_amber_section_i4
1332 :
1333 : ! **************************************************************************************************
1334 : !> \brief Set of Low level subroutines reading section for parmtop
1335 : !> reading 5 arrays of integers of length dim
1336 : !> \param parser ...
1337 : !> \param section ...
1338 : !> \param array1 ...
1339 : !> \param array2 ...
1340 : !> \param array3 ...
1341 : !> \param array4 ...
1342 : !> \param array5 ...
1343 : !> \param dim ...
1344 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1345 : ! **************************************************************************************************
1346 80 : SUBROUTINE rd_amber_section_i5(parser, section, array1, array2, array3, array4, &
1347 80 : array5, dim)
1348 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1349 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1350 : INTEGER, DIMENSION(:) :: array1, array2, array3, array4, array5
1351 : INTEGER, INTENT(IN) :: dim
1352 :
1353 : INTEGER :: i
1354 : LOGICAL :: my_end
1355 :
1356 80 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1357 80 : i = 1
1358 102004 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1359 : !array1
1360 101924 : IF (parser_test_next_token(parser) == "EOL") THEN
1361 152842 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1362 : END IF
1363 101924 : IF (my_end) EXIT
1364 101924 : CALL parser_get_object(parser, array1(i))
1365 : !array2
1366 101924 : IF (parser_test_next_token(parser) == "EOL") THEN
1367 101924 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1368 : END IF
1369 101924 : IF (my_end) EXIT
1370 101924 : CALL parser_get_object(parser, array2(i))
1371 : !array3
1372 101924 : IF (parser_test_next_token(parser) == "EOL") THEN
1373 101924 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1374 : END IF
1375 101924 : IF (my_end) EXIT
1376 101924 : CALL parser_get_object(parser, array3(i))
1377 : !array4
1378 101924 : IF (parser_test_next_token(parser) == "EOL") THEN
1379 101924 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1380 : END IF
1381 101924 : IF (my_end) EXIT
1382 101924 : CALL parser_get_object(parser, array4(i))
1383 : !array5
1384 101924 : IF (parser_test_next_token(parser) == "EOL") THEN
1385 101924 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1386 : END IF
1387 101924 : IF (my_end) EXIT
1388 101924 : CALL parser_get_object(parser, array5(i))
1389 101924 : i = i + 1
1390 : END DO
1391 : ! Trigger end of file aborting
1392 80 : IF (my_end .AND. (i <= dim)) THEN
1393 : CALL cp_abort(__LOCATION__, &
1394 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1395 : END IF
1396 80 : END SUBROUTINE rd_amber_section_i5
1397 :
1398 : ! **************************************************************************************************
1399 : !> \brief Set of Low level subroutines reading section for parmtop
1400 : !> reading 1 array of strings of length dim
1401 : !> \param parser ...
1402 : !> \param section ...
1403 : !> \param array1 ...
1404 : !> \param dim ...
1405 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1406 : ! **************************************************************************************************
1407 48 : SUBROUTINE rd_amber_section_c1(parser, section, array1, dim)
1408 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1409 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1410 : CHARACTER(LEN=default_string_length), DIMENSION(:) :: array1
1411 : INTEGER, INTENT(IN) :: dim
1412 :
1413 : INTEGER :: i
1414 : LOGICAL :: my_end
1415 :
1416 48 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1417 48 : i = 1
1418 101654 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1419 101606 : IF (parser_test_next_token(parser) == "EOL") THEN
1420 106672 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1421 : END IF
1422 101606 : IF (my_end) EXIT
1423 101606 : CALL parser_get_object(parser, array1(i), lower_to_upper=.TRUE.)
1424 101606 : i = i + 1
1425 : END DO
1426 : ! Trigger end of file aborting
1427 48 : IF (my_end .AND. (i <= dim)) THEN
1428 : CALL cp_abort(__LOCATION__, &
1429 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1430 : END IF
1431 48 : END SUBROUTINE rd_amber_section_c1
1432 :
1433 : ! **************************************************************************************************
1434 : !> \brief Set of Low level subroutines reading section for parmtop
1435 : !> reading 1 array of strings of length dim
1436 : !> \param parser ...
1437 : !> \param section ...
1438 : !> \param array1 ...
1439 : !> \param dim ...
1440 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1441 : ! **************************************************************************************************
1442 224 : SUBROUTINE rd_amber_section_r1(parser, section, array1, dim)
1443 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1444 : CHARACTER(LEN=default_string_length), INTENT(IN) :: section
1445 : REAL(KIND=dp), DIMENSION(:) :: array1
1446 : INTEGER, INTENT(IN) :: dim
1447 :
1448 : INTEGER :: i
1449 : LOGICAL :: my_end
1450 :
1451 224 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1452 224 : i = 1
1453 162398 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1454 162174 : IF (parser_test_next_token(parser) == "EOL") THEN
1455 194502 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1456 : END IF
1457 162174 : IF (my_end) EXIT
1458 162174 : CALL parser_get_object(parser, array1(i))
1459 162174 : i = i + 1
1460 : END DO
1461 : ! Trigger end of file aborting
1462 224 : IF (my_end .AND. (i <= dim)) THEN
1463 : CALL cp_abort(__LOCATION__, &
1464 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1465 : END IF
1466 224 : END SUBROUTINE rd_amber_section_r1
1467 :
1468 : ! **************************************************************************************************
1469 : !> \brief Check the version of the AMBER topology file (we can handle from v8 on)
1470 : !> \param parser ...
1471 : !> \param section ...
1472 : !> \param input_format ...
1473 : !> \return ...
1474 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1475 : ! **************************************************************************************************
1476 1560 : FUNCTION get_section_parmtop(parser, section, input_format) RESULT(another_section)
1477 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1478 : CHARACTER(LEN=default_string_length), INTENT(OUT) :: section, input_format
1479 : LOGICAL :: another_section
1480 :
1481 : INTEGER :: end_f, indflag, start_f
1482 : LOGICAL :: found, my_end
1483 :
1484 1560 : CALL parser_search_string(parser, "%FLAG", .TRUE., found, begin_line=.TRUE.)
1485 1560 : IF (found) THEN
1486 : ! section label
1487 1520 : indflag = INDEX(parser%input_line, "%FLAG") + LEN_TRIM("%FLAG")
1488 3040 : DO WHILE (INDEX(parser%input_line(indflag:indflag), " ") /= 0)
1489 1520 : indflag = indflag + 1
1490 : END DO
1491 1520 : section = TRIM(parser%input_line(indflag:))
1492 : ! Input format
1493 1520 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1494 1520 : IF (INDEX(parser%input_line, "%FORMAT") == 0 .OR. my_end) THEN
1495 0 : CPABORT("Expecting %FORMAT. Not found! Abort reading of AMBER topology file!")
1496 : END IF
1497 :
1498 1520 : start_f = INDEX(parser%input_line, "(")
1499 1520 : end_f = INDEX(parser%input_line, ")")
1500 1520 : input_format = parser%input_line(start_f:end_f)
1501 : another_section = .TRUE.
1502 : ELSE
1503 : another_section = .FALSE.
1504 : END IF
1505 1560 : END FUNCTION get_section_parmtop
1506 :
1507 : ! **************************************************************************************************
1508 : !> \brief Check the version of the AMBER topology file (we can handle from v8 on)
1509 : !> \param parser ...
1510 : !> \param output_unit ...
1511 : !> \return ...
1512 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
1513 : ! **************************************************************************************************
1514 40 : FUNCTION check_amber_8_std(parser, output_unit) RESULT(found_AMBER_V8)
1515 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1516 : INTEGER, INTENT(IN) :: output_unit
1517 : LOGICAL :: found_AMBER_V8
1518 :
1519 40 : CALL parser_search_string(parser, "%VERSION ", .TRUE., found_AMBER_V8, begin_line=.TRUE.)
1520 40 : IF (.NOT. found_AMBER_V8) THEN
1521 : CALL cp_abort(__LOCATION__, &
1522 : "This is not an AMBER V.8 PRMTOP format file. Cannot interpret older "// &
1523 0 : "AMBER file formats. ")
1524 : END IF
1525 43 : IF (output_unit > 0) WRITE (output_unit, '(" AMBER_INFO| ",A)') "Amber PrmTop V.8 or greater.", &
1526 6 : TRIM(parser%input_line)
1527 :
1528 40 : END FUNCTION check_amber_8_std
1529 :
1530 : ! **************************************************************************************************
1531 : !> \brief Post processing of forcefield information related to bonds
1532 : !> \param label_a ...
1533 : !> \param label_b ...
1534 : !> \param k ...
1535 : !> \param r0 ...
1536 : !> \param particle_set ...
1537 : !> \param ibond ...
1538 : !> \param nbond ...
1539 : !> \param ib ...
1540 : !> \param jb ...
1541 : !> \param icb ...
1542 : !> \param rk ...
1543 : !> \param req ...
1544 : !> \author Teodoro Laino [tlaino] - 11.2008
1545 : ! **************************************************************************************************
1546 32 : SUBROUTINE post_process_bonds_info(label_a, label_b, k, r0, particle_set, ibond, &
1547 32 : nbond, ib, jb, icb, rk, req)
1548 : CHARACTER(LEN=default_string_length), &
1549 : DIMENSION(:), POINTER :: label_a, label_b
1550 : REAL(KIND=dp), DIMENSION(:), POINTER :: k, r0
1551 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1552 : INTEGER, INTENT(INOUT) :: ibond
1553 : INTEGER, INTENT(IN) :: nbond
1554 : INTEGER, DIMENSION(:), INTENT(IN) :: ib, jb, icb
1555 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rk, req
1556 :
1557 : CHARACTER(len=*), PARAMETER :: routineN = 'post_process_bonds_info'
1558 :
1559 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b
1560 : CHARACTER(LEN=default_string_length), &
1561 32 : ALLOCATABLE, DIMENSION(:, :) :: work_label
1562 : INTEGER :: handle, i
1563 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1564 : LOGICAL :: l_dum
1565 :
1566 32 : CALL timeset(routineN, handle)
1567 32 : IF (nbond /= 0) THEN
1568 90 : ALLOCATE (work_label(2, nbond))
1569 90 : ALLOCATE (iwork(nbond))
1570 56862 : DO i = 1, nbond
1571 56832 : name_atm_a = particle_set(ib(i))%atomic_kind%name
1572 56832 : name_atm_b = particle_set(jb(i))%atomic_kind%name
1573 56832 : l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
1574 56832 : work_label(1, i) = name_atm_a
1575 56862 : work_label(2, i) = name_atm_b
1576 : END DO
1577 30 : CALL sort(work_label, 1, nbond, 1, 2, iwork)
1578 :
1579 30 : ibond = ibond + 1
1580 : ! In case we need more space ... give it up...
1581 30 : IF (ibond > SIZE(label_a)) THEN
1582 4 : CALL reallocate(label_a, 1, INT(buffer_size + ibond*1.5_dp))
1583 4 : CALL reallocate(label_b, 1, INT(buffer_size + ibond*1.5_dp))
1584 4 : CALL reallocate(k, 1, INT(buffer_size + ibond*1.5_dp))
1585 4 : CALL reallocate(r0, 1, INT(buffer_size + ibond*1.5_dp))
1586 : END IF
1587 30 : label_a(ibond) = work_label(1, 1)
1588 30 : label_b(ibond) = work_label(2, 1)
1589 30 : k(ibond) = rk(icb(iwork(1)))
1590 30 : r0(ibond) = req(icb(iwork(1)))
1591 :
1592 56832 : DO i = 2, nbond
1593 56802 : IF ((work_label(1, i) /= label_a(ibond)) .OR. &
1594 30 : (work_label(2, i) /= label_b(ibond))) THEN
1595 1716 : ibond = ibond + 1
1596 : ! In case we need more space ... give it up...
1597 1716 : IF (ibond > SIZE(label_a)) THEN
1598 88 : CALL reallocate(label_a, 1, INT(buffer_size + ibond*1.5_dp))
1599 88 : CALL reallocate(label_b, 1, INT(buffer_size + ibond*1.5_dp))
1600 88 : CALL reallocate(k, 1, INT(buffer_size + ibond*1.5_dp))
1601 88 : CALL reallocate(r0, 1, INT(buffer_size + ibond*1.5_dp))
1602 : END IF
1603 1716 : label_a(ibond) = work_label(1, i)
1604 1716 : label_b(ibond) = work_label(2, i)
1605 1716 : k(ibond) = rk(icb(iwork(i)))
1606 1716 : r0(ibond) = req(icb(iwork(i)))
1607 : END IF
1608 : END DO
1609 :
1610 30 : DEALLOCATE (work_label)
1611 30 : DEALLOCATE (iwork)
1612 : END IF
1613 32 : CALL timestop(handle)
1614 32 : END SUBROUTINE post_process_bonds_info
1615 :
1616 : ! **************************************************************************************************
1617 : !> \brief Post processing of forcefield information related to bends
1618 : !> \param label_a ...
1619 : !> \param label_b ...
1620 : !> \param label_c ...
1621 : !> \param k ...
1622 : !> \param theta0 ...
1623 : !> \param particle_set ...
1624 : !> \param itheta ...
1625 : !> \param ntheta ...
1626 : !> \param it ...
1627 : !> \param jt ...
1628 : !> \param kt ...
1629 : !> \param ict ...
1630 : !> \param tk ...
1631 : !> \param teq ...
1632 : !> \author Teodoro Laino [tlaino] - 11.2008
1633 : ! **************************************************************************************************
1634 32 : SUBROUTINE post_process_bends_info(label_a, label_b, label_c, k, theta0, &
1635 32 : particle_set, itheta, ntheta, it, jt, kt, ict, tk, teq)
1636 : CHARACTER(LEN=default_string_length), &
1637 : DIMENSION(:), POINTER :: label_a, label_b, label_c
1638 : REAL(KIND=dp), DIMENSION(:), POINTER :: k, theta0
1639 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1640 : INTEGER, INTENT(INOUT) :: itheta
1641 : INTEGER, INTENT(IN) :: ntheta
1642 : INTEGER, DIMENSION(:), INTENT(IN) :: it, jt, kt, ict
1643 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tk, teq
1644 :
1645 : CHARACTER(len=*), PARAMETER :: routineN = 'post_process_bends_info'
1646 :
1647 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c
1648 : CHARACTER(LEN=default_string_length), &
1649 32 : ALLOCATABLE, DIMENSION(:, :) :: work_label
1650 : INTEGER :: handle, i
1651 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1652 : LOGICAL :: l_dum
1653 :
1654 32 : CALL timeset(routineN, handle)
1655 32 : IF (ntheta /= 0) THEN
1656 90 : ALLOCATE (work_label(3, ntheta))
1657 90 : ALLOCATE (iwork(ntheta))
1658 45456 : DO i = 1, ntheta
1659 45426 : name_atm_a = particle_set(it(i))%atomic_kind%name
1660 45426 : name_atm_b = particle_set(jt(i))%atomic_kind%name
1661 45426 : name_atm_c = particle_set(kt(i))%atomic_kind%name
1662 : l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, &
1663 45426 : id3=name_atm_c)
1664 45426 : work_label(1, i) = name_atm_a
1665 45426 : work_label(2, i) = name_atm_b
1666 45456 : work_label(3, i) = name_atm_c
1667 : END DO
1668 :
1669 30 : CALL sort(work_label, 1, ntheta, 1, 3, iwork)
1670 :
1671 30 : itheta = itheta + 1
1672 : ! In case we need more space ... give it up...
1673 30 : IF (itheta > SIZE(label_a)) THEN
1674 2 : CALL reallocate(label_a, 1, INT(buffer_size + itheta*1.5_dp))
1675 2 : CALL reallocate(label_b, 1, INT(buffer_size + itheta*1.5_dp))
1676 2 : CALL reallocate(label_c, 1, INT(buffer_size + itheta*1.5_dp))
1677 2 : CALL reallocate(k, 1, INT(buffer_size + itheta*1.5_dp))
1678 2 : CALL reallocate(theta0, 1, INT(buffer_size + itheta*1.5_dp))
1679 : END IF
1680 30 : label_a(itheta) = work_label(1, 1)
1681 30 : label_b(itheta) = work_label(2, 1)
1682 30 : label_c(itheta) = work_label(3, 1)
1683 30 : k(itheta) = tk(ict(iwork(1)))
1684 30 : theta0(itheta) = teq(ict(iwork(1)))
1685 :
1686 45426 : DO i = 2, ntheta
1687 : IF ((work_label(1, i) /= label_a(itheta)) .OR. &
1688 45396 : (work_label(2, i) /= label_b(itheta)) .OR. &
1689 30 : (work_label(3, i) /= label_c(itheta))) THEN
1690 3642 : itheta = itheta + 1
1691 : ! In case we need more space ... give it up...
1692 3642 : IF (itheta > SIZE(label_a)) THEN
1693 110 : CALL reallocate(label_a, 1, INT(buffer_size + itheta*1.5_dp))
1694 110 : CALL reallocate(label_b, 1, INT(buffer_size + itheta*1.5_dp))
1695 110 : CALL reallocate(label_c, 1, INT(buffer_size + itheta*1.5_dp))
1696 110 : CALL reallocate(k, 1, INT(buffer_size + itheta*1.5_dp))
1697 110 : CALL reallocate(theta0, 1, INT(buffer_size + itheta*1.5_dp))
1698 : END IF
1699 3642 : label_a(itheta) = work_label(1, i)
1700 3642 : label_b(itheta) = work_label(2, i)
1701 3642 : label_c(itheta) = work_label(3, i)
1702 3642 : k(itheta) = tk(ict(iwork(i)))
1703 3642 : theta0(itheta) = teq(ict(iwork(i)))
1704 : END IF
1705 : END DO
1706 :
1707 30 : DEALLOCATE (work_label)
1708 30 : DEALLOCATE (iwork)
1709 : END IF
1710 32 : CALL timestop(handle)
1711 32 : END SUBROUTINE post_process_bends_info
1712 :
1713 : ! **************************************************************************************************
1714 : !> \brief Post processing of forcefield information related to torsions
1715 : !> \param label_a ...
1716 : !> \param label_b ...
1717 : !> \param label_c ...
1718 : !> \param label_d ...
1719 : !> \param k ...
1720 : !> \param m ...
1721 : !> \param phi0 ...
1722 : !> \param particle_set ...
1723 : !> \param iphi ...
1724 : !> \param nphi ...
1725 : !> \param ip ...
1726 : !> \param jp ...
1727 : !> \param kp ...
1728 : !> \param lp ...
1729 : !> \param icp ...
1730 : !> \param pk ...
1731 : !> \param pn ...
1732 : !> \param phase ...
1733 : !> \author Teodoro Laino [tlaino] - 11.2008
1734 : ! **************************************************************************************************
1735 32 : SUBROUTINE post_process_torsions_info(label_a, label_b, label_c, label_d, k, &
1736 32 : m, phi0, particle_set, iphi, nphi, ip, jp, kp, lp, icp, pk, pn, phase)
1737 : CHARACTER(LEN=default_string_length), &
1738 : DIMENSION(:), POINTER :: label_a, label_b, label_c, label_d
1739 : REAL(KIND=dp), DIMENSION(:), POINTER :: k
1740 : INTEGER, DIMENSION(:), POINTER :: m
1741 : REAL(KIND=dp), DIMENSION(:), POINTER :: phi0
1742 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1743 : INTEGER, INTENT(INOUT) :: iphi
1744 : INTEGER, INTENT(IN) :: nphi
1745 : INTEGER, DIMENSION(:), INTENT(IN) :: ip, jp, kp, lp, icp
1746 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: pk, pn, phase
1747 :
1748 : CHARACTER(len=*), PARAMETER :: routineN = 'post_process_torsions_info'
1749 :
1750 : CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c, &
1751 : name_atm_d
1752 : CHARACTER(LEN=default_string_length), &
1753 32 : ALLOCATABLE, DIMENSION(:, :) :: work_label
1754 : INTEGER :: handle, i
1755 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1756 : LOGICAL :: l_dum
1757 :
1758 32 : CALL timeset(routineN, handle)
1759 32 : IF (nphi /= 0) THEN
1760 84 : ALLOCATE (work_label(6, nphi))
1761 84 : ALLOCATE (iwork(nphi))
1762 50488 : DO i = 1, nphi
1763 50460 : name_atm_a = particle_set(ip(i))%atomic_kind%name
1764 50460 : name_atm_b = particle_set(jp(i))%atomic_kind%name
1765 50460 : name_atm_c = particle_set(kp(i))%atomic_kind%name
1766 50460 : name_atm_d = particle_set(lp(i))%atomic_kind%name
1767 : l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, &
1768 50460 : id3=name_atm_c, id4=name_atm_d)
1769 50460 : work_label(1, i) = name_atm_a
1770 50460 : work_label(2, i) = name_atm_b
1771 50460 : work_label(3, i) = name_atm_c
1772 50460 : work_label(4, i) = name_atm_d
1773 : ! Phase and multiplicity must be kept into account
1774 : ! for the ordering of the torsions
1775 50460 : work_label(5, i) = TRIM(ADJUSTL(cp_to_string(phase(icp(i)))))
1776 50488 : work_label(6, i) = TRIM(ADJUSTL(cp_to_string(pn(icp(i)))))
1777 : END DO
1778 :
1779 28 : CALL sort(work_label, 1, nphi, 1, 6, iwork)
1780 :
1781 28 : iphi = iphi + 1
1782 : ! In case we need more space ... give it up...
1783 28 : IF (iphi > SIZE(label_a)) THEN
1784 0 : CALL reallocate(label_a, 1, INT(buffer_size + iphi*1.5_dp))
1785 0 : CALL reallocate(label_b, 1, INT(buffer_size + iphi*1.5_dp))
1786 0 : CALL reallocate(label_c, 1, INT(buffer_size + iphi*1.5_dp))
1787 0 : CALL reallocate(label_d, 1, INT(buffer_size + iphi*1.5_dp))
1788 0 : CALL reallocate(k, 1, INT(buffer_size + iphi*1.5_dp))
1789 0 : CALL reallocate(m, 1, INT(buffer_size + iphi*1.5_dp))
1790 0 : CALL reallocate(phi0, 1, INT(buffer_size + iphi*1.5_dp))
1791 : END IF
1792 28 : label_a(iphi) = work_label(1, 1)
1793 28 : label_b(iphi) = work_label(2, 1)
1794 28 : label_c(iphi) = work_label(3, 1)
1795 28 : label_d(iphi) = work_label(4, 1)
1796 28 : k(iphi) = pk(icp(iwork(1)))
1797 28 : m(iphi) = NINT(pn(icp(iwork(1))))
1798 28 : IF (m(iphi) - pn(icp(iwork(1))) > EPSILON(1.0_dp)) THEN
1799 0 : CPABORT("Non-integer torsions not supported")
1800 : END IF
1801 :
1802 28 : phi0(iphi) = phase(icp(iwork(1)))
1803 :
1804 50460 : DO i = 2, nphi
1805 : ! We don't consider the possibility that a torsion can have same
1806 : ! phase, periodicity but different value of k.. in this case the
1807 : ! potential should be summed-up
1808 : IF ((work_label(1, i) /= label_a(iphi)) .OR. &
1809 : (work_label(2, i) /= label_b(iphi)) .OR. &
1810 : (work_label(3, i) /= label_c(iphi)) .OR. &
1811 : (work_label(4, i) /= label_d(iphi)) .OR. &
1812 50432 : (pn(icp(iwork(i))) /= m(iphi)) .OR. &
1813 28 : (phase(icp(iwork(i))) /= phi0(iphi))) THEN
1814 10096 : iphi = iphi + 1
1815 : ! In case we need more space ... give it up...
1816 10096 : IF (iphi > SIZE(label_a)) THEN
1817 138 : CALL reallocate(label_a, 1, INT(buffer_size + iphi*1.5_dp))
1818 138 : CALL reallocate(label_b, 1, INT(buffer_size + iphi*1.5_dp))
1819 138 : CALL reallocate(label_c, 1, INT(buffer_size + iphi*1.5_dp))
1820 138 : CALL reallocate(label_d, 1, INT(buffer_size + iphi*1.5_dp))
1821 138 : CALL reallocate(k, 1, INT(buffer_size + iphi*1.5_dp))
1822 138 : CALL reallocate(m, 1, INT(buffer_size + iphi*1.5_dp))
1823 138 : CALL reallocate(phi0, 1, INT(buffer_size + iphi*1.5_dp))
1824 : END IF
1825 10096 : label_a(iphi) = work_label(1, i)
1826 10096 : label_b(iphi) = work_label(2, i)
1827 10096 : label_c(iphi) = work_label(3, i)
1828 10096 : label_d(iphi) = work_label(4, i)
1829 10096 : k(iphi) = pk(icp(iwork(i)))
1830 10096 : m(iphi) = NINT(pn(icp(iwork(i))))
1831 10096 : IF (m(iphi) - pn(icp(iwork(i))) > EPSILON(1.0_dp)) THEN
1832 0 : CPABORT("Non-integer torsions not supported")
1833 : END IF
1834 10096 : phi0(iphi) = phase(icp(iwork(i)))
1835 : END IF
1836 : END DO
1837 :
1838 28 : DEALLOCATE (work_label)
1839 28 : DEALLOCATE (iwork)
1840 : END IF
1841 32 : CALL timestop(handle)
1842 32 : END SUBROUTINE post_process_torsions_info
1843 :
1844 : ! **************************************************************************************************
1845 : !> \brief Post processing of forcefield information related to Lennard-Jones
1846 : !> \param atom_label ...
1847 : !> \param eps ...
1848 : !> \param sigma ...
1849 : !> \param particle_set ...
1850 : !> \param ntypes ...
1851 : !> \param nsize ...
1852 : !> \param iac ...
1853 : !> \param ico ...
1854 : !> \param cn1 ...
1855 : !> \param cn2 ...
1856 : !> \param natom ...
1857 : !> \author Teodoro Laino [tlaino] - 11.2008
1858 : ! **************************************************************************************************
1859 16 : SUBROUTINE post_process_LJ_info(atom_label, eps, sigma, particle_set, &
1860 16 : ntypes, nsize, iac, ico, cn1, cn2, natom)
1861 : CHARACTER(LEN=default_string_length), &
1862 : DIMENSION(:), POINTER :: atom_label
1863 : REAL(KIND=dp), DIMENSION(:), POINTER :: eps, sigma
1864 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1865 : INTEGER, INTENT(IN) :: ntypes
1866 : INTEGER, INTENT(INOUT) :: nsize
1867 : INTEGER, DIMENSION(:), INTENT(IN) :: iac, ico
1868 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cn1, cn2
1869 : INTEGER, INTENT(IN) :: natom
1870 :
1871 : CHARACTER(len=*), PARAMETER :: routineN = 'post_process_LJ_info'
1872 :
1873 : CHARACTER(LEN=default_string_length) :: name_atm_a
1874 : CHARACTER(LEN=default_string_length), &
1875 16 : ALLOCATABLE, DIMENSION(:) :: work_label
1876 : INTEGER :: handle, i
1877 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1878 : LOGICAL :: check, l_dum
1879 : REAL(KIND=dp) :: F12, F6, my_eps, my_sigma, sigma6
1880 :
1881 16 : CALL timeset(routineN, handle)
1882 48 : ALLOCATE (work_label(natom))
1883 48 : ALLOCATE (iwork(natom))
1884 78544 : DO i = 1, natom
1885 78528 : name_atm_a = particle_set(i)%atomic_kind%name
1886 78528 : l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a)
1887 78544 : work_label(i) = name_atm_a
1888 : END DO
1889 16 : CALL sort(work_label, natom, iwork)
1890 :
1891 16 : nsize = nsize + 1
1892 16 : IF (nsize > SIZE(atom_label)) THEN
1893 0 : CALL reallocate(atom_label, 1, INT(buffer_size + nsize*1.5_dp))
1894 0 : CALL reallocate(eps, 1, INT(buffer_size + nsize*1.5_dp))
1895 0 : CALL reallocate(sigma, 1, INT(buffer_size + nsize*1.5_dp))
1896 : END IF
1897 16 : F12 = cn1(ico(ntypes*(iac(iwork(1)) - 1) + iac(iwork(1))))
1898 16 : F6 = cn2(ico(ntypes*(iac(iwork(1)) - 1) + iac(iwork(1))))
1899 16 : check = (F6 == 0.0_dp) .EQV. (F12 == 0.0_dp)
1900 16 : CPASSERT(check)
1901 16 : my_sigma = 0.0_dp
1902 16 : my_eps = 0.0_dp
1903 16 : IF (F6 /= 0.0_dp) THEN
1904 16 : sigma6 = (2.0_dp*F12/F6)
1905 16 : my_sigma = sigma6**(1.0_dp/6.0_dp)
1906 16 : my_eps = F6/(2.0_dp*sigma6)
1907 : END IF
1908 16 : atom_label(nsize) = work_label(1)
1909 16 : sigma(nsize) = my_sigma/2.0_dp
1910 16 : eps(nsize) = my_eps
1911 :
1912 78528 : DO i = 2, natom
1913 78528 : IF (work_label(i) == atom_label(nsize)) THEN
1914 77020 : IF (iac(iwork(i)) /= iac(iwork(i - 1))) THEN
1915 0 : CPABORT("AMBER kind collision: "//TRIM(work_label(i)))
1916 : END IF
1917 : ELSE
1918 1492 : nsize = nsize + 1
1919 : ! In case we need more space ... give it up...
1920 1492 : IF (nsize > SIZE(atom_label)) THEN
1921 90 : CALL reallocate(atom_label, 1, INT(buffer_size + nsize*1.5_dp))
1922 90 : CALL reallocate(eps, 1, INT(buffer_size + nsize*1.5_dp))
1923 90 : CALL reallocate(sigma, 1, INT(buffer_size + nsize*1.5_dp))
1924 : END IF
1925 1492 : F12 = cn1(ico(ntypes*(iac(iwork(i)) - 1) + iac(iwork(i))))
1926 1492 : F6 = cn2(ico(ntypes*(iac(iwork(i)) - 1) + iac(iwork(i))))
1927 1492 : check = (F6 == 0.0_dp) .EQV. (F12 == 0.0_dp)
1928 1492 : CPASSERT(check)
1929 1492 : my_sigma = 0.0_dp
1930 1492 : my_eps = 0.0_dp
1931 1492 : IF (F6 /= 0.0_dp) THEN
1932 1468 : sigma6 = (2.0_dp*F12/F6)
1933 1468 : my_sigma = sigma6**(1.0_dp/6.0_dp)
1934 1468 : my_eps = F6/(2.0_dp*sigma6)
1935 : END IF
1936 1492 : atom_label(nsize) = work_label(i)
1937 1492 : sigma(nsize) = my_sigma/2.0_dp
1938 1492 : eps(nsize) = my_eps
1939 : END IF
1940 : END DO
1941 :
1942 16 : DEALLOCATE (work_label)
1943 16 : DEALLOCATE (iwork)
1944 16 : CALL timestop(handle)
1945 16 : END SUBROUTINE post_process_LJ_info
1946 :
1947 : END MODULE topology_amber
1948 :
|
