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 22 : 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 22 : NULLIFY (logger)
300 22 : CALL timeset(routineN, handle)
301 22 : logger => cp_get_default_logger()
302 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/AMBER_INFO", &
303 22 : extension=".subsysLog")
304 :
305 22 : atom_info => topology%atom_info
306 22 : 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 22 : atom_info=atom_info, conn_info=conn_info)
311 :
312 : ! Molnames have been internally generated
313 22 : topology%molname_generated = .TRUE.
314 :
315 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
316 22 : "PRINT%TOPOLOGY_INFO/AMBER_INFO")
317 22 : CALL timestop(handle)
318 22 : 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 36 : 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 72 : 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 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iac, ib, ibh, icb, icbh, ico, icp, icph, &
391 36 : ict, icth, ip, iph, ipres, it, ith, &
392 36 : iwork, jb, jbh, jp, jph, jt, jth, kp, &
393 36 : kph, kt, kth, lp, lph
394 36 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: full_torsions
395 : LOGICAL :: check, valid_format
396 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: asol, bsol, cn1, cn2, phase, pk, pn, &
397 36 : req, rk, teq, tk
398 : TYPE(cp_parser_type) :: parser
399 :
400 36 : CALL timeset(routineN, handle)
401 36 : IF (output_unit > 0) WRITE (output_unit, '(/,A)') " AMBER_INFO| Reading Amber Topology File: "// &
402 3 : TRIM(filename)
403 36 : CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.TRUE.)
404 36 : valid_format = check_amber_8_std(parser, output_unit)
405 36 : IF (valid_format) THEN
406 1412 : DO WHILE (get_section_parmtop(parser, section, input_format))
407 36 : SELECT CASE (TRIM(section))
408 : CASE ("TITLE")
409 : ! Who cares about the title?
410 36 : CYCLE
411 : CASE ("POINTERS")
412 36 : 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 36 : natom = info(1)
416 36 : ntypes = info(2)
417 36 : nbonh = info(3)
418 36 : mbona = info(4)
419 36 : ntheth = info(5)
420 36 : mtheta = info(6)
421 36 : nphih = info(7)
422 36 : mphia = info(8)
423 36 : nhparm = info(9)
424 36 : nparm = info(10)
425 36 : nnb = info(11)
426 36 : nres = info(12)
427 36 : nbona = info(13)
428 36 : ntheta = info(14)
429 36 : nphia = info(15)
430 36 : numbnd = info(16)
431 36 : numang = info(17)
432 36 : nptra = info(18)
433 36 : natyp = info(19)
434 36 : nphb = info(20)
435 36 : ifpert = info(21)
436 36 : nbper = info(22)
437 36 : ngper = info(23)
438 36 : ndper = info(24)
439 36 : mbper = info(25)
440 36 : mgper = info(26)
441 36 : mdper = info(27)
442 36 : ifbox = info(28)
443 36 : nmxrs = info(29)
444 36 : ifcap = info(30)
445 36 : numextra = info(31)
446 :
447 : ! Print some info if requested
448 36 : 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 36 : IF (do_connectivity) THEN
459 22 : check = PRESENT(atom_info) .AND. PRESENT(conn_info)
460 22 : CPASSERT(check)
461 22 : natom_prev = 0
462 22 : IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname)
463 : ! Allocate for extracting connectivity infos
464 66 : ALLOCATE (labres(nres))
465 66 : ALLOCATE (ipres(nres))
466 : END IF
467 36 : IF (do_forcefield) THEN
468 : ! Allocate for extracting forcefield infos
469 42 : ALLOCATE (iac(natom))
470 42 : ALLOCATE (ico(ntypes*ntypes))
471 42 : ALLOCATE (rk(numbnd))
472 28 : ALLOCATE (req(numbnd))
473 42 : ALLOCATE (tk(numang))
474 28 : ALLOCATE (teq(numang))
475 40 : ALLOCATE (pk(nptra))
476 26 : ALLOCATE (pn(nptra))
477 26 : ALLOCATE (phase(nptra))
478 42 : ALLOCATE (cn1(ntypes*(ntypes + 1)/2))
479 28 : ALLOCATE (cn2(ntypes*(ntypes + 1)/2))
480 28 : ALLOCATE (asol(ntypes*(ntypes + 1)/2))
481 28 : ALLOCATE (bsol(ntypes*(ntypes + 1)/2))
482 : END IF
483 : ! Always Allocate
484 108 : ALLOCATE (ibh(nbonh))
485 72 : ALLOCATE (jbh(nbonh))
486 72 : ALLOCATE (icbh(nbonh))
487 104 : ALLOCATE (ib(nbona))
488 68 : ALLOCATE (jb(nbona))
489 68 : ALLOCATE (icb(nbona))
490 108 : ALLOCATE (ith(ntheth))
491 72 : ALLOCATE (jth(ntheth))
492 72 : ALLOCATE (kth(ntheth))
493 72 : ALLOCATE (icth(ntheth))
494 104 : ALLOCATE (it(ntheta))
495 68 : ALLOCATE (jt(ntheta))
496 68 : ALLOCATE (kt(ntheta))
497 68 : ALLOCATE (ict(ntheta))
498 104 : ALLOCATE (iph(nphih))
499 68 : ALLOCATE (jph(nphih))
500 68 : ALLOCATE (kph(nphih))
501 68 : ALLOCATE (lph(nphih))
502 68 : ALLOCATE (icph(nphih))
503 104 : ALLOCATE (ip(nphia))
504 68 : ALLOCATE (jp(nphia))
505 68 : ALLOCATE (kp(nphia))
506 68 : ALLOCATE (lp(nphia))
507 68 : ALLOCATE (icp(nphia))
508 : CASE ("ATOM_NAME")
509 : ! Atom names are just ignored according the CP2K philosophy
510 36 : CYCLE
511 : CASE ("AMBER_ATOM_TYPE")
512 36 : IF (.NOT. do_connectivity) CYCLE
513 22 : CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom)
514 66 : ALLOCATE (strtmp_a(natom))
515 22 : CALL rd_amber_section(parser, section, strtmp_a, natom)
516 78860 : DO i = 1, natom
517 78860 : atom_info%id_atmname(natom_prev + i) = str2id(strtmp_a(i))
518 : END DO
519 22 : DEALLOCATE (strtmp_a)
520 : CASE ("CHARGE")
521 36 : IF (.NOT. do_connectivity) CYCLE
522 22 : CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom)
523 22 : CALL rd_amber_section(parser, section, atom_info%atm_charge(natom_prev + 1:), natom)
524 : ! Convert charges into atomic units
525 78860 : atom_info%atm_charge(natom_prev + 1:) = atom_info%atm_charge(natom_prev + 1:)/amber_conv_charge
526 : CASE ("MASS")
527 36 : IF (.NOT. do_connectivity) CYCLE
528 22 : CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom)
529 22 : CALL rd_amber_section(parser, section, atom_info%atm_mass(natom_prev + 1:), natom)
530 : CASE ("RESIDUE_LABEL")
531 36 : IF (.NOT. do_connectivity) CYCLE
532 22 : CALL reallocate(atom_info%id_resname, 1, natom_prev + natom)
533 22 : CALL rd_amber_section(parser, section, labres, nres)
534 : CASE ("RESIDUE_POINTER")
535 36 : IF (.NOT. do_connectivity) CYCLE
536 22 : CALL reallocate(atom_info%resid, 1, natom_prev + natom)
537 22 : CALL rd_amber_section(parser, section, ipres, nres)
538 : CASE ("ATOM_TYPE_INDEX")
539 36 : IF (.NOT. do_forcefield) CYCLE
540 14 : CALL rd_amber_section(parser, section, iac, natom)
541 : CASE ("NONBONDED_PARM_INDEX")
542 36 : IF (.NOT. do_forcefield) CYCLE
543 14 : CALL rd_amber_section(parser, section, ico, ntypes**2)
544 : CASE ("BOND_FORCE_CONSTANT")
545 36 : IF (.NOT. do_forcefield) CYCLE
546 14 : CALL rd_amber_section(parser, section, rk, numbnd)
547 : CASE ("BOND_EQUIL_VALUE")
548 36 : IF (.NOT. do_forcefield) CYCLE
549 14 : CALL rd_amber_section(parser, section, req, numbnd)
550 : CASE ("ANGLE_FORCE_CONSTANT")
551 36 : IF (.NOT. do_forcefield) CYCLE
552 14 : CALL rd_amber_section(parser, section, tk, numang)
553 : CASE ("ANGLE_EQUIL_VALUE")
554 36 : IF (.NOT. do_forcefield) CYCLE
555 14 : CALL rd_amber_section(parser, section, teq, numang)
556 : CASE ("DIHEDRAL_FORCE_CONSTANT")
557 36 : IF (.NOT. do_forcefield) CYCLE
558 14 : CALL rd_amber_section(parser, section, pk, nptra)
559 14 : IF (nptra <= 0) CYCLE
560 : ! Save raw values
561 12 : IF (ASSOCIATED(amb_info%raw_torsion_k)) DEALLOCATE (amb_info%raw_torsion_k)
562 358 : ALLOCATE (amb_info%raw_torsion_k(nptra), source=pk)
563 : CASE ("DIHEDRAL_PERIODICITY")
564 36 : IF (.NOT. do_forcefield) CYCLE
565 14 : CALL rd_amber_section(parser, section, pn, nptra)
566 14 : IF (nptra <= 0) CYCLE
567 : ! Save raw values
568 12 : IF (ASSOCIATED(amb_info%raw_torsion_m)) DEALLOCATE (amb_info%raw_torsion_m)
569 358 : ALLOCATE (amb_info%raw_torsion_m(nptra), source=pn)
570 : CASE ("DIHEDRAL_PHASE")
571 36 : IF (.NOT. do_forcefield) CYCLE
572 14 : CALL rd_amber_section(parser, section, phase, nptra)
573 14 : IF (nptra <= 0) CYCLE
574 : ! Save raw values
575 12 : IF (ASSOCIATED(amb_info%raw_torsion_phi0)) DEALLOCATE (amb_info%raw_torsion_phi0)
576 358 : ALLOCATE (amb_info%raw_torsion_phi0(nptra), source=phase)
577 : CASE ("LENNARD_JONES_ACOEF")
578 36 : IF (.NOT. do_forcefield) CYCLE
579 14 : CALL rd_amber_section(parser, section, cn1, ntypes*(ntypes + 1)/2)
580 : CASE ("LENNARD_JONES_BCOEF")
581 36 : IF (.NOT. do_forcefield) CYCLE
582 14 : CALL rd_amber_section(parser, section, cn2, ntypes*(ntypes + 1)/2)
583 : CASE ("HBOND_ACOEF")
584 36 : IF (.NOT. do_forcefield) CYCLE
585 14 : CALL rd_amber_section(parser, section, asol, nphb)
586 : CASE ("HBOND_BCOEF")
587 36 : IF (.NOT. do_forcefield) CYCLE
588 14 : 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 36 : CALL rd_amber_section(parser, section, ibh, jbh, icbh, nbonh)
592 : ! Conver to an atomic index
593 100028 : ibh(:) = ibh(:)/3 + 1
594 100028 : jbh(:) = jbh(:)/3 + 1
595 : CASE ("BONDS_WITHOUT_HYDROGEN")
596 : ! We always need to parse this information both for connectivity and forcefields
597 36 : CALL rd_amber_section(parser, section, ib, jb, icb, nbona)
598 : ! Conver to an atomic index
599 14022 : ib(:) = ib(:)/3 + 1
600 14022 : jb(:) = jb(:)/3 + 1
601 : CASE ("ANGLES_INC_HYDROGEN")
602 : ! We always need to parse this information both for connectivity and forcefields
603 36 : CALL rd_amber_section(parser, section, ith, jth, kth, icth, ntheth)
604 : ! Conver to an atomic index
605 72486 : ith(:) = ith(:)/3 + 1
606 72486 : jth(:) = jth(:)/3 + 1
607 72486 : kth(:) = kth(:)/3 + 1
608 : CASE ("ANGLES_WITHOUT_HYDROGEN")
609 : ! We always need to parse this information both for connectivity and forcefields
610 36 : CALL rd_amber_section(parser, section, it, jt, kt, ict, ntheta)
611 : ! Conver to an atomic index
612 18954 : it(:) = it(:)/3 + 1
613 18954 : jt(:) = jt(:)/3 + 1
614 18954 : kt(:) = kt(:)/3 + 1
615 : CASE ("DIHEDRALS_INC_HYDROGEN")
616 : ! We always need to parse this information both for connectivity and forcefields
617 36 : CALL rd_amber_section(parser, section, iph, jph, kph, lph, icph, nphih)
618 : ! Conver to an atomic index
619 56580 : iph(:) = iph(:)/3 + 1
620 56580 : jph(:) = jph(:)/3 + 1
621 56580 : kph(:) = ABS(kph(:))/3 + 1
622 56580 : lph(:) = ABS(lph(:))/3 + 1
623 : CASE ("DIHEDRALS_WITHOUT_HYDROGEN")
624 : ! We always need to parse this information both for connectivity and forcefields
625 36 : CALL rd_amber_section(parser, section, ip, jp, kp, lp, icp, nphia)
626 : ! Conver to an atomic index
627 45272 : ip(:) = ip(:)/3 + 1
628 45272 : jp(:) = jp(:)/3 + 1
629 45272 : kp(:) = ABS(kp(:))/3 + 1
630 46662 : 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 36 : IF (do_forcefield .AND. (nphih + nphia > 0)) THEN
637 12 : IF (ASSOCIATED(amb_info%raw_torsion_id)) DEALLOCATE (amb_info%raw_torsion_id)
638 36 : ALLOCATE (amb_info%raw_torsion_id(5, nphih + nphia))
639 28078 : DO i = 1, nphih
640 28066 : amb_info%raw_torsion_id(1, i) = iph(i)
641 28066 : amb_info%raw_torsion_id(2, i) = jph(i)
642 28066 : amb_info%raw_torsion_id(3, i) = kph(i)
643 28066 : amb_info%raw_torsion_id(4, i) = lph(i)
644 28078 : amb_info%raw_torsion_id(5, i) = icph(i)
645 : END DO
646 22334 : DO i = 1, nphia
647 22322 : amb_info%raw_torsion_id(1, nphih + i) = ip(i)
648 22322 : amb_info%raw_torsion_id(2, nphih + i) = jp(i)
649 22322 : amb_info%raw_torsion_id(3, nphih + i) = kp(i)
650 22322 : amb_info%raw_torsion_id(4, nphih + i) = lp(i)
651 22334 : 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 36 : IF (do_connectivity) THEN
658 22 : CALL timeset(TRIM(routineN)//"_connectivity", handle2)
659 : ! ----------------------------------------------------------
660 : ! Conform Amber Names with CHARMM convention (kind<->charge)
661 : ! ----------------------------------------------------------
662 66 : ALLOCATE (isymbl(natom))
663 66 : ALLOCATE (iwork(natom))
664 :
665 78860 : DO i = 1, SIZE(isymbl)
666 78860 : isymbl(i) = id2str(atom_info%id_atmname(natom_prev + i))
667 : END DO
668 :
669 : ! Sort atom names + charges and identify unique types
670 22 : CALL sort(isymbl, natom, iwork)
671 :
672 22 : istart = 1
673 78838 : DO i = 2, natom
674 78838 : IF (TRIM(isymbl(i)) /= TRIM(isymbl(istart))) THEN
675 228 : CALL conform_atom_type_low(isymbl, iwork, i, istart, atom_info%atm_charge(natom_prev + 1:))
676 228 : istart = i
677 : END IF
678 : END DO
679 22 : 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 78860 : DO i = 1, natom
683 78860 : 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 22730 : DO i = 1, nres - 1
690 123776 : atom_info%id_resname(natom_prev + ipres(i):natom_prev + ipres(i + 1)) = str2id(s2s(labres(i)))
691 123798 : atom_info%resid(natom_prev + ipres(i):natom_prev + ipres(i + 1)) = i
692 : END DO
693 500 : atom_info%id_resname(natom_prev + ipres(i):natom_prev + natom) = str2id(s2s(labres(i)))
694 500 : atom_info%resid(natom_prev + ipres(i):natom_prev + natom) = i
695 :
696 : ! Deallocate when extracting connectivity infos
697 22 : DEALLOCATE (iwork)
698 22 : DEALLOCATE (isymbl)
699 22 : DEALLOCATE (labres)
700 22 : DEALLOCATE (ipres)
701 :
702 : ! ----------------------------------------------------------
703 : ! Copy connectivity
704 : ! ----------------------------------------------------------
705 : ! BONDS
706 22 : nbond_prev = 0
707 22 : IF (ASSOCIATED(conn_info%bond_a)) nbond_prev = SIZE(conn_info%bond_a)
708 :
709 22 : CALL reallocate(conn_info%bond_a, 1, nbond_prev + nbonh + nbona)
710 22 : CALL reallocate(conn_info%bond_b, 1, nbond_prev + nbonh + nbona)
711 50078 : DO i = 1, nbonh
712 50056 : index_now = nbond_prev + i
713 50056 : conn_info%bond_a(index_now) = natom_prev + ibh(i)
714 50078 : conn_info%bond_b(index_now) = natom_prev + jbh(i)
715 : END DO
716 7144 : DO i = 1, nbona
717 7122 : index_now = nbond_prev + i + nbonh
718 7122 : conn_info%bond_a(index_now) = natom_prev + ib(i)
719 7144 : conn_info%bond_b(index_now) = natom_prev + jb(i)
720 : END DO
721 :
722 : ! ANGLES
723 22 : ntheta_prev = 0
724 22 : IF (ASSOCIATED(conn_info%theta_a)) ntheta_prev = SIZE(conn_info%theta_a)
725 :
726 22 : CALL reallocate(conn_info%theta_a, 1, ntheta_prev + ntheth + ntheta)
727 22 : CALL reallocate(conn_info%theta_b, 1, ntheta_prev + ntheth + ntheta)
728 22 : CALL reallocate(conn_info%theta_c, 1, ntheta_prev + ntheth + ntheta)
729 36368 : DO i = 1, ntheth
730 36346 : index_now = ntheta_prev + i
731 36346 : conn_info%theta_a(index_now) = natom_prev + ith(i)
732 36346 : conn_info%theta_b(index_now) = natom_prev + jth(i)
733 36368 : conn_info%theta_c(index_now) = natom_prev + kth(i)
734 : END DO
735 9672 : DO i = 1, ntheta
736 9650 : index_now = ntheta_prev + i + ntheth
737 9650 : conn_info%theta_a(index_now) = natom_prev + it(i)
738 9650 : conn_info%theta_b(index_now) = natom_prev + jt(i)
739 9672 : 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 22 : nphi_prev = 0
746 22 : IF (ASSOCIATED(conn_info%phi_a)) nphi_prev = SIZE(conn_info%phi_a)
747 :
748 22 : CALL reallocate(conn_info%phi_a, 1, nphi_prev + nphih + nphia)
749 22 : CALL reallocate(conn_info%phi_b, 1, nphi_prev + nphih + nphia)
750 22 : CALL reallocate(conn_info%phi_c, 1, nphi_prev + nphih + nphia)
751 22 : CALL reallocate(conn_info%phi_d, 1, nphi_prev + nphih + nphia)
752 :
753 22 : IF (nphih + nphia /= 0) THEN
754 60 : ALLOCATE (full_torsions(4, nphih + nphia))
755 60 : ALLOCATE (iwork(nphih + nphia))
756 :
757 28498 : DO i = 1, nphih
758 28478 : full_torsions(1, i) = iph(i)
759 28478 : full_torsions(2, i) = jph(i)
760 28478 : full_torsions(3, i) = kph(i)
761 28498 : full_torsions(4, i) = lph(i)
762 : END DO
763 22934 : DO i = 1, nphia
764 22914 : full_torsions(1, nphih + i) = ip(i)
765 22914 : full_torsions(2, nphih + i) = jp(i)
766 22914 : full_torsions(3, nphih + i) = kp(i)
767 22934 : full_torsions(4, nphih + i) = lp(i)
768 : END DO
769 20 : CALL sort(full_torsions, 1, nphih + nphia, 1, 4, iwork)
770 :
771 20 : unique_torsions = nphi_prev + 1
772 20 : conn_info%phi_a(unique_torsions) = natom_prev + full_torsions(1, 1)
773 20 : conn_info%phi_b(unique_torsions) = natom_prev + full_torsions(2, 1)
774 20 : conn_info%phi_c(unique_torsions) = natom_prev + full_torsions(3, 1)
775 20 : conn_info%phi_d(unique_torsions) = natom_prev + full_torsions(4, 1)
776 51392 : 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 51372 : (full_torsions(3, i) /= full_torsions(3, i - 1)) .OR. &
780 20 : (full_torsions(4, i) /= full_torsions(4, i - 1))) THEN
781 37586 : unique_torsions = unique_torsions + 1
782 37586 : conn_info%phi_a(unique_torsions) = natom_prev + full_torsions(1, i)
783 37586 : conn_info%phi_b(unique_torsions) = natom_prev + full_torsions(2, i)
784 37586 : conn_info%phi_c(unique_torsions) = natom_prev + full_torsions(3, i)
785 37586 : conn_info%phi_d(unique_torsions) = natom_prev + full_torsions(4, i)
786 : END IF
787 : END DO
788 20 : CALL reallocate(conn_info%phi_a, 1, unique_torsions)
789 20 : CALL reallocate(conn_info%phi_b, 1, unique_torsions)
790 20 : CALL reallocate(conn_info%phi_c, 1, unique_torsions)
791 20 : CALL reallocate(conn_info%phi_d, 1, unique_torsions)
792 :
793 20 : DEALLOCATE (full_torsions)
794 20 : DEALLOCATE (iwork)
795 : END IF
796 : ! IMPROPERS
797 22 : CALL reallocate(conn_info%impr_a, 1, 0)
798 22 : CALL reallocate(conn_info%impr_b, 1, 0)
799 22 : CALL reallocate(conn_info%impr_c, 1, 0)
800 22 : CALL reallocate(conn_info%impr_d, 1, 0)
801 :
802 : ! ----------------------------------------------------------
803 : ! Generate molecule names
804 : ! ----------------------------------------------------------
805 22 : CALL reallocate(atom_info%id_molname, 1, natom_prev + natom)
806 78860 : 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 22 : atom_info%id_molname(natom_prev + 1:natom_prev + natom))
809 44 : CALL timestop(handle2)
810 : END IF
811 :
812 : ! Extracts force fields info from the AMBER topology file
813 36 : IF (do_forcefield) THEN
814 14 : CALL timeset(TRIM(routineN)//"_forcefield", handle2)
815 : ! ----------------------------------------------------------
816 : ! Force Fields informations related to bonds
817 : ! ----------------------------------------------------------
818 14 : CALL reallocate(amb_info%bond_a, 1, buffer_size)
819 14 : CALL reallocate(amb_info%bond_b, 1, buffer_size)
820 14 : CALL reallocate(amb_info%bond_k, 1, buffer_size)
821 14 : CALL reallocate(amb_info%bond_r0, 1, buffer_size)
822 14 : 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 14 : 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 14 : nbona, ib, jb, icb, rk, req)
831 : ! Shrink arrays size to the minimal request
832 14 : CALL reallocate(amb_info%bond_a, 1, nsize)
833 14 : CALL reallocate(amb_info%bond_b, 1, nsize)
834 14 : CALL reallocate(amb_info%bond_k, 1, nsize)
835 14 : CALL reallocate(amb_info%bond_r0, 1, nsize)
836 :
837 : ! ----------------------------------------------------------
838 : ! Force Fields informations related to bends
839 : ! ----------------------------------------------------------
840 14 : CALL reallocate(amb_info%bend_a, 1, buffer_size)
841 14 : CALL reallocate(amb_info%bend_b, 1, buffer_size)
842 14 : CALL reallocate(amb_info%bend_c, 1, buffer_size)
843 14 : CALL reallocate(amb_info%bend_k, 1, buffer_size)
844 14 : CALL reallocate(amb_info%bend_theta0, 1, buffer_size)
845 14 : 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 14 : 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 14 : particle_set, nsize, ntheta, it, jt, kt, ict, tk, teq)
854 : ! Shrink arrays size to the minimal request
855 14 : CALL reallocate(amb_info%bend_a, 1, nsize)
856 14 : CALL reallocate(amb_info%bend_b, 1, nsize)
857 14 : CALL reallocate(amb_info%bend_c, 1, nsize)
858 14 : CALL reallocate(amb_info%bend_k, 1, nsize)
859 14 : 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 14 : CALL reallocate(amb_info%torsion_a, 1, buffer_size)
867 14 : CALL reallocate(amb_info%torsion_b, 1, buffer_size)
868 14 : CALL reallocate(amb_info%torsion_c, 1, buffer_size)
869 14 : CALL reallocate(amb_info%torsion_d, 1, buffer_size)
870 14 : CALL reallocate(amb_info%torsion_k, 1, buffer_size)
871 14 : CALL reallocate(amb_info%torsion_m, 1, buffer_size)
872 14 : CALL reallocate(amb_info%torsion_phi0, 1, buffer_size)
873 14 : 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 14 : 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 14 : nphia, ip, jp, kp, lp, icp, pk, pn, phase)
884 : ! Shrink arrays size to the minimal request
885 14 : CALL reallocate(amb_info%torsion_a, 1, nsize)
886 14 : CALL reallocate(amb_info%torsion_b, 1, nsize)
887 14 : CALL reallocate(amb_info%torsion_c, 1, nsize)
888 14 : CALL reallocate(amb_info%torsion_d, 1, nsize)
889 14 : CALL reallocate(amb_info%torsion_k, 1, nsize)
890 14 : CALL reallocate(amb_info%torsion_m, 1, nsize)
891 14 : CALL reallocate(amb_info%torsion_phi0, 1, nsize)
892 :
893 : ! Sort dihedral metadata for faster lookup
894 14 : IF (nphih + nphia /= 0) THEN
895 36 : ALLOCATE (iwork(nphih + nphia))
896 12 : CALL sort(amb_info%raw_torsion_id, 1, nphih + nphia, 1, 5, iwork)
897 12 : DEALLOCATE (iwork)
898 : END IF
899 :
900 : ! ----------------------------------------------------------
901 : ! Post process of LJ parameters
902 : ! ----------------------------------------------------------
903 14 : CALL reallocate(amb_info%nonbond_a, 1, buffer_size)
904 14 : CALL reallocate(amb_info%nonbond_eps, 1, buffer_size)
905 14 : CALL reallocate(amb_info%nonbond_rmin2, 1, buffer_size)
906 :
907 14 : 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 14 : cn1, cn2, natom)
911 :
912 : ! Shrink arrays size to the minimal request
913 14 : CALL reallocate(amb_info%nonbond_a, 1, nsize)
914 14 : CALL reallocate(amb_info%nonbond_eps, 1, nsize)
915 14 : CALL reallocate(amb_info%nonbond_rmin2, 1, nsize)
916 :
917 : ! Deallocate at the end of the dirty job
918 14 : DEALLOCATE (iac)
919 14 : DEALLOCATE (ico)
920 14 : DEALLOCATE (rk)
921 14 : DEALLOCATE (req)
922 14 : DEALLOCATE (tk)
923 14 : DEALLOCATE (teq)
924 14 : DEALLOCATE (pk)
925 14 : DEALLOCATE (pn)
926 14 : DEALLOCATE (phase)
927 14 : DEALLOCATE (cn1)
928 14 : DEALLOCATE (cn2)
929 14 : DEALLOCATE (asol)
930 14 : DEALLOCATE (bsol)
931 14 : CALL timestop(handle2)
932 : END IF
933 : ! Always Deallocate
934 36 : DEALLOCATE (ibh)
935 36 : DEALLOCATE (jbh)
936 36 : DEALLOCATE (icbh)
937 36 : DEALLOCATE (ib)
938 36 : DEALLOCATE (jb)
939 36 : DEALLOCATE (icb)
940 36 : DEALLOCATE (ith)
941 36 : DEALLOCATE (jth)
942 36 : DEALLOCATE (kth)
943 36 : DEALLOCATE (icth)
944 36 : DEALLOCATE (it)
945 36 : DEALLOCATE (jt)
946 36 : DEALLOCATE (kt)
947 36 : DEALLOCATE (ict)
948 36 : DEALLOCATE (iph)
949 36 : DEALLOCATE (jph)
950 36 : DEALLOCATE (kph)
951 36 : DEALLOCATE (lph)
952 36 : DEALLOCATE (icph)
953 36 : DEALLOCATE (ip)
954 36 : DEALLOCATE (jp)
955 36 : DEALLOCATE (kp)
956 36 : DEALLOCATE (lp)
957 36 : DEALLOCATE (icp)
958 36 : CALL parser_release(parser)
959 36 : CALL timestop(handle)
960 36 : RETURN
961 : ! Output info Format
962 : 1000 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 144 : 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 250 : 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 250 : INTEGER, DIMENSION(:), POINTER :: cindx, lindx
992 : REAL(KIND=dp) :: ctmp
993 250 : REAL(KIND=dp), DIMENSION(:), POINTER :: cwork
994 :
995 250 : CALL timeset(routineN, handle)
996 250 : iend = i - 1
997 250 : isize = iend - istart + 1
998 750 : ALLOCATE (cwork(isize))
999 750 : ALLOCATE (lindx(isize))
1000 500 : ALLOCATE (cindx(isize))
1001 250 : ind = 0
1002 79088 : DO k = istart, iend
1003 78838 : ind = ind + 1
1004 78838 : cwork(ind) = charges(iwork(k))
1005 79088 : lindx(ind) = k
1006 : END DO
1007 250 : CALL sort(cwork, isize, cindx)
1008 :
1009 250 : ctmp = cwork(1)
1010 250 : counter = 1
1011 78838 : DO k = 2, isize
1012 78838 : IF (cwork(k) /= ctmp) THEN
1013 1408 : counter = counter + 1
1014 1408 : ctmp = cwork(k)
1015 : END IF
1016 : END DO
1017 250 : IF (counter /= 1) THEN
1018 148 : counter = 1
1019 148 : kstart = 1
1020 148 : ctmp = cwork(1)
1021 12762 : DO k = 2, isize
1022 12762 : IF (cwork(k) /= ctmp) THEN
1023 : kend = k - 1
1024 12348 : DO j = kstart, kend
1025 10940 : gind = lindx(cindx(j))
1026 12348 : isymbl(gind) = TRIM(isymbl(gind))//ADJUSTL(cp_to_string(counter))
1027 : END DO
1028 1408 : counter = counter + 1
1029 1408 : ctmp = cwork(k)
1030 1408 : kstart = k
1031 : END IF
1032 : END DO
1033 : kend = k - 1
1034 1970 : DO j = kstart, kend
1035 1822 : gind = lindx(cindx(j))
1036 1970 : isymbl(gind) = TRIM(isymbl(gind))//ADJUSTL(cp_to_string(counter))
1037 : END DO
1038 : END IF
1039 250 : DEALLOCATE (cwork)
1040 250 : DEALLOCATE (lindx)
1041 250 : DEALLOCATE (cindx)
1042 250 : CALL timestop(handle)
1043 250 : 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 86 : 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 86 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1064 86 : i = 1
1065 104356 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1066 104270 : IF (parser_test_next_token(parser) == "EOL") THEN
1067 114666 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1068 : END IF
1069 104270 : IF (my_end) EXIT
1070 104270 : CALL parser_get_object(parser, array1(i))
1071 104270 : i = i + 1
1072 : END DO
1073 : ! Trigger end of file aborting
1074 86 : IF (my_end .AND. (i <= dim)) THEN
1075 : CALL cp_abort(__LOCATION__, &
1076 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1077 : END IF
1078 86 : 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 72 : 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 72 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1101 72 : i = 1
1102 114050 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1103 : !array1
1104 113978 : IF (parser_test_next_token(parser) == "EOL") THEN
1105 125336 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1106 : END IF
1107 113978 : IF (my_end) EXIT
1108 113978 : CALL parser_get_object(parser, array1(i))
1109 : !array2
1110 113978 : IF (parser_test_next_token(parser) == "EOL") THEN
1111 125398 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1112 : END IF
1113 113978 : IF (my_end) EXIT
1114 113978 : CALL parser_get_object(parser, array2(i))
1115 : !array3
1116 113978 : IF (parser_test_next_token(parser) == "EOL") THEN
1117 125356 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1118 : END IF
1119 113978 : IF (my_end) EXIT
1120 113978 : CALL parser_get_object(parser, array3(i))
1121 113978 : i = i + 1
1122 : END DO
1123 : ! Trigger end of file aborting
1124 72 : IF (my_end .AND. (i <= dim)) THEN
1125 : CALL cp_abort(__LOCATION__, &
1126 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1127 : END IF
1128 72 : 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 72 : 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 72 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1152 72 : i = 1
1153 91440 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1154 : !array1
1155 91368 : IF (parser_test_next_token(parser) == "EOL") THEN
1156 109606 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1157 : END IF
1158 91368 : IF (my_end) EXIT
1159 91368 : CALL parser_get_object(parser, array1(i))
1160 : !array2
1161 91368 : IF (parser_test_next_token(parser) == "EOL") THEN
1162 91368 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1163 : END IF
1164 91368 : IF (my_end) EXIT
1165 91368 : CALL parser_get_object(parser, array2(i))
1166 : !array3
1167 91368 : IF (parser_test_next_token(parser) == "EOL") THEN
1168 109634 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1169 : END IF
1170 91368 : IF (my_end) EXIT
1171 91368 : CALL parser_get_object(parser, array3(i))
1172 : !array4
1173 91368 : IF (parser_test_next_token(parser) == "EOL") THEN
1174 91368 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1175 : END IF
1176 91368 : IF (my_end) EXIT
1177 91368 : CALL parser_get_object(parser, array4(i))
1178 91368 : i = i + 1
1179 : END DO
1180 : ! Trigger end of file aborting
1181 72 : IF (my_end .AND. (i <= dim)) THEN
1182 : CALL cp_abort(__LOCATION__, &
1183 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1184 : END IF
1185 72 : 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 72 : SUBROUTINE rd_amber_section_i5(parser, section, array1, array2, array3, array4, &
1201 72 : 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 72 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1211 72 : i = 1
1212 101852 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1213 : !array1
1214 101780 : IF (parser_test_next_token(parser) == "EOL") THEN
1215 152634 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1216 : END IF
1217 101780 : IF (my_end) EXIT
1218 101780 : CALL parser_get_object(parser, array1(i))
1219 : !array2
1220 101780 : IF (parser_test_next_token(parser) == "EOL") THEN
1221 101780 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1222 : END IF
1223 101780 : IF (my_end) EXIT
1224 101780 : CALL parser_get_object(parser, array2(i))
1225 : !array3
1226 101780 : IF (parser_test_next_token(parser) == "EOL") THEN
1227 101780 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1228 : END IF
1229 101780 : IF (my_end) EXIT
1230 101780 : CALL parser_get_object(parser, array3(i))
1231 : !array4
1232 101780 : IF (parser_test_next_token(parser) == "EOL") THEN
1233 101780 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1234 : END IF
1235 101780 : IF (my_end) EXIT
1236 101780 : CALL parser_get_object(parser, array4(i))
1237 : !array5
1238 101780 : IF (parser_test_next_token(parser) == "EOL") THEN
1239 101780 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1240 : END IF
1241 101780 : IF (my_end) EXIT
1242 101780 : CALL parser_get_object(parser, array5(i))
1243 101780 : i = i + 1
1244 : END DO
1245 : ! Trigger end of file aborting
1246 72 : IF (my_end .AND. (i <= dim)) THEN
1247 : CALL cp_abort(__LOCATION__, &
1248 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1249 : END IF
1250 72 : 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 44 : 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 44 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1271 44 : i = 1
1272 101612 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1273 101568 : IF (parser_test_next_token(parser) == "EOL") THEN
1274 106634 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1275 : END IF
1276 101568 : IF (my_end) EXIT
1277 101568 : CALL parser_get_object(parser, array1(i), lower_to_upper=.TRUE.)
1278 101568 : i = i + 1
1279 : END DO
1280 : ! Trigger end of file aborting
1281 44 : IF (my_end .AND. (i <= dim)) THEN
1282 : CALL cp_abort(__LOCATION__, &
1283 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1284 : END IF
1285 44 : 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 198 : 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 198 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1306 198 : i = 1
1307 162032 : DO WHILE ((i <= dim) .AND. (.NOT. my_end))
1308 161834 : IF (parser_test_next_token(parser) == "EOL") THEN
1309 194108 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1310 : END IF
1311 161834 : IF (my_end) EXIT
1312 161834 : CALL parser_get_object(parser, array1(i))
1313 161834 : i = i + 1
1314 : END DO
1315 : ! Trigger end of file aborting
1316 198 : IF (my_end .AND. (i <= dim)) THEN
1317 : CALL cp_abort(__LOCATION__, &
1318 0 : "End of file while reading section "//TRIM(section)//" in amber topology file!")
1319 : END IF
1320 198 : 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 1412 : 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 1412 : CALL parser_search_string(parser, "%FLAG", .TRUE., found, begin_line=.TRUE.)
1339 1412 : IF (found) THEN
1340 : ! section label
1341 1376 : indflag = INDEX(parser%input_line, "%FLAG") + LEN_TRIM("%FLAG")
1342 2752 : DO WHILE (INDEX(parser%input_line(indflag:indflag), " ") /= 0)
1343 1376 : indflag = indflag + 1
1344 : END DO
1345 1376 : section = TRIM(parser%input_line(indflag:))
1346 : ! Input format
1347 1376 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1348 1376 : IF (INDEX(parser%input_line, "%FORMAT") == 0 .OR. my_end) THEN
1349 0 : CPABORT("Expecting %FORMAT. Not found! Abort reading of AMBER topology file!")
1350 : END IF
1351 :
1352 1376 : start_f = INDEX(parser%input_line, "(")
1353 1376 : end_f = INDEX(parser%input_line, ")")
1354 1376 : input_format = parser%input_line(start_f:end_f)
1355 : another_section = .TRUE.
1356 : ELSE
1357 : another_section = .FALSE.
1358 : END IF
1359 1412 : 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 36 : 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 36 : CALL parser_search_string(parser, "%VERSION ", .TRUE., found_AMBER_V8, begin_line=.TRUE.)
1374 36 : 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 0 : "AMBER file formats. ")
1378 : END IF
1379 39 : IF (output_unit > 0) WRITE (output_unit, '(" AMBER_INFO| ",A)') "Amber PrmTop V.8 or greater.", &
1380 6 : TRIM(parser%input_line)
1381 :
1382 36 : 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 28 : SUBROUTINE post_process_bonds_info(label_a, label_b, k, r0, particle_set, ibond, &
1401 28 : 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 28 : ALLOCATABLE, DIMENSION(:, :) :: work_label
1416 : INTEGER :: handle, i
1417 28 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1418 : LOGICAL :: l_dum
1419 :
1420 28 : CALL timeset(routineN, handle)
1421 28 : IF (nbond /= 0) THEN
1422 78 : ALLOCATE (work_label(2, nbond))
1423 78 : ALLOCATE (iwork(nbond))
1424 56826 : DO i = 1, nbond
1425 56800 : name_atm_a = particle_set(ib(i))%atomic_kind%name
1426 56800 : name_atm_b = particle_set(jb(i))%atomic_kind%name
1427 56800 : l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b)
1428 56800 : work_label(1, i) = name_atm_a
1429 56826 : work_label(2, i) = name_atm_b
1430 : END DO
1431 26 : CALL sort(work_label, 1, nbond, 1, 2, iwork)
1432 :
1433 26 : ibond = ibond + 1
1434 : ! In case we need more space ... give it up...
1435 26 : IF (ibond > SIZE(label_a)) THEN
1436 2 : CALL reallocate(label_a, 1, INT(buffer_size + ibond*1.5_dp))
1437 2 : CALL reallocate(label_b, 1, INT(buffer_size + ibond*1.5_dp))
1438 2 : CALL reallocate(k, 1, INT(buffer_size + ibond*1.5_dp))
1439 2 : CALL reallocate(r0, 1, INT(buffer_size + ibond*1.5_dp))
1440 : END IF
1441 26 : label_a(ibond) = work_label(1, 1)
1442 26 : label_b(ibond) = work_label(2, 1)
1443 26 : k(ibond) = rk(icb(iwork(1)))
1444 26 : r0(ibond) = req(icb(iwork(1)))
1445 :
1446 56800 : DO i = 2, nbond
1447 56774 : IF ((work_label(1, i) /= label_a(ibond)) .OR. &
1448 26 : (work_label(2, i) /= label_b(ibond))) THEN
1449 1698 : ibond = ibond + 1
1450 : ! In case we need more space ... give it up...
1451 1698 : IF (ibond > SIZE(label_a)) THEN
1452 84 : CALL reallocate(label_a, 1, INT(buffer_size + ibond*1.5_dp))
1453 84 : CALL reallocate(label_b, 1, INT(buffer_size + ibond*1.5_dp))
1454 84 : CALL reallocate(k, 1, INT(buffer_size + ibond*1.5_dp))
1455 84 : CALL reallocate(r0, 1, INT(buffer_size + ibond*1.5_dp))
1456 : END IF
1457 1698 : label_a(ibond) = work_label(1, i)
1458 1698 : label_b(ibond) = work_label(2, i)
1459 1698 : k(ibond) = rk(icb(iwork(i)))
1460 1698 : r0(ibond) = req(icb(iwork(i)))
1461 : END IF
1462 : END DO
1463 :
1464 26 : DEALLOCATE (work_label)
1465 26 : DEALLOCATE (iwork)
1466 : END IF
1467 28 : CALL timestop(handle)
1468 28 : 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 28 : SUBROUTINE post_process_bends_info(label_a, label_b, label_c, k, theta0, &
1489 28 : 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 28 : ALLOCATABLE, DIMENSION(:, :) :: work_label
1504 : INTEGER :: handle, i
1505 28 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1506 : LOGICAL :: l_dum
1507 :
1508 28 : CALL timeset(routineN, handle)
1509 28 : IF (ntheta /= 0) THEN
1510 78 : ALLOCATE (work_label(3, ntheta))
1511 78 : ALLOCATE (iwork(ntheta))
1512 45398 : DO i = 1, ntheta
1513 45372 : name_atm_a = particle_set(it(i))%atomic_kind%name
1514 45372 : name_atm_b = particle_set(jt(i))%atomic_kind%name
1515 45372 : 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 45372 : id3=name_atm_c)
1518 45372 : work_label(1, i) = name_atm_a
1519 45372 : work_label(2, i) = name_atm_b
1520 45398 : work_label(3, i) = name_atm_c
1521 : END DO
1522 :
1523 26 : CALL sort(work_label, 1, ntheta, 1, 3, iwork)
1524 :
1525 26 : itheta = itheta + 1
1526 : ! In case we need more space ... give it up...
1527 26 : IF (itheta > SIZE(label_a)) THEN
1528 2 : CALL reallocate(label_a, 1, INT(buffer_size + itheta*1.5_dp))
1529 2 : CALL reallocate(label_b, 1, INT(buffer_size + itheta*1.5_dp))
1530 2 : CALL reallocate(label_c, 1, INT(buffer_size + itheta*1.5_dp))
1531 2 : CALL reallocate(k, 1, INT(buffer_size + itheta*1.5_dp))
1532 2 : CALL reallocate(theta0, 1, INT(buffer_size + itheta*1.5_dp))
1533 : END IF
1534 26 : label_a(itheta) = work_label(1, 1)
1535 26 : label_b(itheta) = work_label(2, 1)
1536 26 : label_c(itheta) = work_label(3, 1)
1537 26 : k(itheta) = tk(ict(iwork(1)))
1538 26 : theta0(itheta) = teq(ict(iwork(1)))
1539 :
1540 45372 : DO i = 2, ntheta
1541 : IF ((work_label(1, i) /= label_a(itheta)) .OR. &
1542 45346 : (work_label(2, i) /= label_b(itheta)) .OR. &
1543 26 : (work_label(3, i) /= label_c(itheta))) THEN
1544 3610 : itheta = itheta + 1
1545 : ! In case we need more space ... give it up...
1546 3610 : IF (itheta > SIZE(label_a)) THEN
1547 102 : CALL reallocate(label_a, 1, INT(buffer_size + itheta*1.5_dp))
1548 102 : CALL reallocate(label_b, 1, INT(buffer_size + itheta*1.5_dp))
1549 102 : CALL reallocate(label_c, 1, INT(buffer_size + itheta*1.5_dp))
1550 102 : CALL reallocate(k, 1, INT(buffer_size + itheta*1.5_dp))
1551 102 : CALL reallocate(theta0, 1, INT(buffer_size + itheta*1.5_dp))
1552 : END IF
1553 3610 : label_a(itheta) = work_label(1, i)
1554 3610 : label_b(itheta) = work_label(2, i)
1555 3610 : label_c(itheta) = work_label(3, i)
1556 3610 : k(itheta) = tk(ict(iwork(i)))
1557 3610 : theta0(itheta) = teq(ict(iwork(i)))
1558 : END IF
1559 : END DO
1560 :
1561 26 : DEALLOCATE (work_label)
1562 26 : DEALLOCATE (iwork)
1563 : END IF
1564 28 : CALL timestop(handle)
1565 28 : 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 28 : SUBROUTINE post_process_torsions_info(label_a, label_b, label_c, label_d, k, &
1590 28 : 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 28 : ALLOCATABLE, DIMENSION(:, :) :: work_label
1608 : INTEGER :: handle, i
1609 28 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1610 : LOGICAL :: l_dum
1611 :
1612 28 : CALL timeset(routineN, handle)
1613 28 : IF (nphi /= 0) THEN
1614 72 : ALLOCATE (work_label(6, nphi))
1615 72 : ALLOCATE (iwork(nphi))
1616 50412 : DO i = 1, nphi
1617 50388 : name_atm_a = particle_set(ip(i))%atomic_kind%name
1618 50388 : name_atm_b = particle_set(jp(i))%atomic_kind%name
1619 50388 : name_atm_c = particle_set(kp(i))%atomic_kind%name
1620 50388 : 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 50388 : id3=name_atm_c, id4=name_atm_d)
1623 50388 : work_label(1, i) = name_atm_a
1624 50388 : work_label(2, i) = name_atm_b
1625 50388 : work_label(3, i) = name_atm_c
1626 50388 : work_label(4, i) = name_atm_d
1627 : ! Phase and multiplicity must be kept into account
1628 : ! for the ordering of the torsions
1629 50388 : work_label(5, i) = TRIM(ADJUSTL(cp_to_string(phase(icp(i)))))
1630 50412 : work_label(6, i) = TRIM(ADJUSTL(cp_to_string(pn(icp(i)))))
1631 : END DO
1632 :
1633 24 : CALL sort(work_label, 1, nphi, 1, 6, iwork)
1634 :
1635 24 : iphi = iphi + 1
1636 : ! In case we need more space ... give it up...
1637 24 : IF (iphi > SIZE(label_a)) THEN
1638 0 : CALL reallocate(label_a, 1, INT(buffer_size + iphi*1.5_dp))
1639 0 : CALL reallocate(label_b, 1, INT(buffer_size + iphi*1.5_dp))
1640 0 : CALL reallocate(label_c, 1, INT(buffer_size + iphi*1.5_dp))
1641 0 : CALL reallocate(label_d, 1, INT(buffer_size + iphi*1.5_dp))
1642 0 : CALL reallocate(k, 1, INT(buffer_size + iphi*1.5_dp))
1643 0 : CALL reallocate(m, 1, INT(buffer_size + iphi*1.5_dp))
1644 0 : CALL reallocate(phi0, 1, INT(buffer_size + iphi*1.5_dp))
1645 : END IF
1646 24 : label_a(iphi) = work_label(1, 1)
1647 24 : label_b(iphi) = work_label(2, 1)
1648 24 : label_c(iphi) = work_label(3, 1)
1649 24 : label_d(iphi) = work_label(4, 1)
1650 24 : k(iphi) = pk(icp(iwork(1)))
1651 24 : m(iphi) = NINT(pn(icp(iwork(1))))
1652 24 : IF (m(iphi) - pn(icp(iwork(1))) > EPSILON(1.0_dp)) THEN
1653 0 : CPABORT("Non-integer torsions not supported")
1654 : END IF
1655 :
1656 24 : phi0(iphi) = phase(icp(iwork(1)))
1657 :
1658 50388 : 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 50364 : (pn(icp(iwork(i))) /= m(iphi)) .OR. &
1667 24 : (phase(icp(iwork(i))) /= phi0(iphi))) THEN
1668 10058 : iphi = iphi + 1
1669 : ! In case we need more space ... give it up...
1670 10058 : IF (iphi > SIZE(label_a)) THEN
1671 130 : CALL reallocate(label_a, 1, INT(buffer_size + iphi*1.5_dp))
1672 130 : CALL reallocate(label_b, 1, INT(buffer_size + iphi*1.5_dp))
1673 130 : CALL reallocate(label_c, 1, INT(buffer_size + iphi*1.5_dp))
1674 130 : CALL reallocate(label_d, 1, INT(buffer_size + iphi*1.5_dp))
1675 130 : CALL reallocate(k, 1, INT(buffer_size + iphi*1.5_dp))
1676 130 : CALL reallocate(m, 1, INT(buffer_size + iphi*1.5_dp))
1677 130 : CALL reallocate(phi0, 1, INT(buffer_size + iphi*1.5_dp))
1678 : END IF
1679 10058 : label_a(iphi) = work_label(1, i)
1680 10058 : label_b(iphi) = work_label(2, i)
1681 10058 : label_c(iphi) = work_label(3, i)
1682 10058 : label_d(iphi) = work_label(4, i)
1683 10058 : k(iphi) = pk(icp(iwork(i)))
1684 10058 : m(iphi) = NINT(pn(icp(iwork(i))))
1685 10058 : IF (m(iphi) - pn(icp(iwork(i))) > EPSILON(1.0_dp)) THEN
1686 0 : CPABORT("Non-integer torsions not supported")
1687 : END IF
1688 10058 : phi0(iphi) = phase(icp(iwork(i)))
1689 : END IF
1690 : END DO
1691 :
1692 24 : DEALLOCATE (work_label)
1693 24 : DEALLOCATE (iwork)
1694 : END IF
1695 28 : CALL timestop(handle)
1696 28 : 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 14 : SUBROUTINE post_process_LJ_info(atom_label, eps, sigma, particle_set, &
1714 14 : 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 14 : ALLOCATABLE, DIMENSION(:) :: work_label
1730 : INTEGER :: handle, i
1731 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1732 : LOGICAL :: check, l_dum
1733 : REAL(KIND=dp) :: F12, F6, my_eps, my_sigma, sigma6
1734 :
1735 14 : CALL timeset(routineN, handle)
1736 42 : ALLOCATE (work_label(natom))
1737 42 : ALLOCATE (iwork(natom))
1738 78508 : DO i = 1, natom
1739 78494 : name_atm_a = particle_set(i)%atomic_kind%name
1740 78494 : l_dum = qmmm_ff_precond_only_qm(id1=name_atm_a)
1741 78508 : work_label(i) = name_atm_a
1742 : END DO
1743 14 : CALL sort(work_label, natom, iwork)
1744 :
1745 14 : nsize = nsize + 1
1746 14 : IF (nsize > SIZE(atom_label)) THEN
1747 0 : CALL reallocate(atom_label, 1, INT(buffer_size + nsize*1.5_dp))
1748 0 : CALL reallocate(eps, 1, INT(buffer_size + nsize*1.5_dp))
1749 0 : CALL reallocate(sigma, 1, INT(buffer_size + nsize*1.5_dp))
1750 : END IF
1751 14 : F12 = cn1(ico(ntypes*(iac(iwork(1)) - 1) + iac(iwork(1))))
1752 14 : F6 = cn2(ico(ntypes*(iac(iwork(1)) - 1) + iac(iwork(1))))
1753 14 : check = (F6 == 0.0_dp) .EQV. (F12 == 0.0_dp)
1754 14 : CPASSERT(check)
1755 14 : my_sigma = 0.0_dp
1756 14 : my_eps = 0.0_dp
1757 14 : IF (F6 /= 0.0_dp) THEN
1758 14 : sigma6 = (2.0_dp*F12/F6)
1759 14 : my_sigma = sigma6**(1.0_dp/6.0_dp)
1760 14 : my_eps = F6/(2.0_dp*sigma6)
1761 : END IF
1762 14 : atom_label(nsize) = work_label(1)
1763 14 : sigma(nsize) = my_sigma/2.0_dp
1764 14 : eps(nsize) = my_eps
1765 :
1766 78494 : DO i = 2, natom
1767 78494 : IF (work_label(i) /= atom_label(nsize)) THEN
1768 1446 : nsize = nsize + 1
1769 : ! In case we need more space ... give it up...
1770 1446 : IF (nsize > SIZE(atom_label)) THEN
1771 84 : CALL reallocate(atom_label, 1, INT(buffer_size + nsize*1.5_dp))
1772 84 : CALL reallocate(eps, 1, INT(buffer_size + nsize*1.5_dp))
1773 84 : CALL reallocate(sigma, 1, INT(buffer_size + nsize*1.5_dp))
1774 : END IF
1775 1446 : F12 = cn1(ico(ntypes*(iac(iwork(i)) - 1) + iac(iwork(i))))
1776 1446 : F6 = cn2(ico(ntypes*(iac(iwork(i)) - 1) + iac(iwork(i))))
1777 1446 : check = (F6 == 0.0_dp) .EQV. (F12 == 0.0_dp)
1778 1446 : CPASSERT(check)
1779 1446 : my_sigma = 0.0_dp
1780 1446 : my_eps = 0.0_dp
1781 1446 : IF (F6 /= 0.0_dp) THEN
1782 1422 : sigma6 = (2.0_dp*F12/F6)
1783 1422 : my_sigma = sigma6**(1.0_dp/6.0_dp)
1784 1422 : my_eps = F6/(2.0_dp*sigma6)
1785 : END IF
1786 1446 : atom_label(nsize) = work_label(i)
1787 1446 : sigma(nsize) = my_sigma/2.0_dp
1788 1446 : eps(nsize) = my_eps
1789 : END IF
1790 : END DO
1791 :
1792 14 : DEALLOCATE (work_label)
1793 14 : DEALLOCATE (iwork)
1794 14 : CALL timestop(handle)
1795 14 : END SUBROUTINE post_process_LJ_info
1796 :
1797 : END MODULE topology_amber
1798 :
|
