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