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