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 : !> \par History
10 : !> Subroutine input_torsions changed (DG) 05-Dec-2000
11 : !> Output formats changed (DG) 05-Dec-2000
12 : !> JGH (26-01-2002) : force field parameters stored in tables, not in
13 : !> matrices. Input changed to have parameters labeled by the position
14 : !> and not atom pairs (triples etc)
15 : !> Teo (11.2005) : Moved all information on force field pair_potential to
16 : !> a much lighter memory structure
17 : !> Teo 09.2006 : Split all routines force_field I/O in a separate file
18 : !> 10.2008 Teodoro Laino [tlaino] : splitted from force_fields_input in order
19 : !> to collect in a unique module only the functions to read external FF
20 : !> \author CJM
21 : ! **************************************************************************************************
22 : MODULE force_fields_ext
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_type
25 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
26 : cp_print_key_unit_nr
27 : USE cp_parser_methods, ONLY: parser_get_next_line,&
28 : parser_get_object,&
29 : parser_search_string,&
30 : parser_test_next_token
31 : USE cp_parser_types, ONLY: cp_parser_type,&
32 : parser_create,&
33 : parser_release
34 : USE cp_units, ONLY: cp_unit_to_cp2k
35 : USE force_field_kind_types, ONLY: do_ff_g96
36 : USE force_field_types, ONLY: amber_info_type,&
37 : charmm_info_type,&
38 : force_field_type,&
39 : gromos_info_type
40 : USE input_section_types, ONLY: section_vals_type
41 : USE kinds, ONLY: default_string_length,&
42 : dp
43 : USE memory_utilities, ONLY: reallocate
44 : USE message_passing, ONLY: mp_para_env_type
45 : USE particle_types, ONLY: particle_type
46 : USE string_utilities, ONLY: uppercase
47 : USE topology_amber, ONLY: rdparm_amber_8
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_ext'
53 :
54 : PRIVATE
55 : PUBLIC :: read_force_field_charmm, &
56 : read_force_field_amber, &
57 : read_force_field_gromos
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Reads the GROMOS force_field
63 : !> \param ff_type ...
64 : !> \param para_env ...
65 : !> \param mm_section ...
66 : !> \author ikuo
67 : ! **************************************************************************************************
68 72 : SUBROUTINE read_force_field_gromos(ff_type, para_env, mm_section)
69 :
70 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
71 : TYPE(mp_para_env_type), POINTER :: para_env
72 : TYPE(section_vals_type), POINTER :: mm_section
73 :
74 : CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_gromos'
75 : REAL(KIND=dp), PARAMETER :: ekt = 2.5_dp
76 :
77 : CHARACTER(LEN=default_string_length) :: label
78 : CHARACTER(LEN=default_string_length), &
79 : DIMENSION(21) :: avail_section
80 12 : CHARACTER(LEN=default_string_length), POINTER :: namearray(:)
81 : INTEGER :: handle, iatom, icon, itemp, itype, iw, &
82 : jatom, ncon, ntype, offset
83 : LOGICAL :: found
84 : REAL(KIND=dp) :: cosphi0, cost2, csq, sdet
85 : TYPE(cp_logger_type), POINTER :: logger
86 : TYPE(cp_parser_type) :: parser
87 : TYPE(gromos_info_type), POINTER :: gro_info
88 :
89 12 : CALL timeset(routineN, handle)
90 12 : NULLIFY (logger)
91 12 : logger => cp_get_default_logger()
92 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
93 12 : extension=".mmLog")
94 :
95 12 : avail_section(1) = "TITLE"
96 12 : avail_section(2) = "TOPPHYSCON"
97 12 : avail_section(3) = "TOPVERSION"
98 12 : avail_section(4) = "ATOMTYPENAME"
99 12 : avail_section(5) = "RESNAME"
100 12 : avail_section(6) = "SOLUTEATOM"
101 12 : avail_section(7) = "BONDTYPE"
102 12 : avail_section(8) = "BONDH"
103 12 : avail_section(9) = "BOND"
104 12 : avail_section(10) = "BONDANGLETYPE"
105 12 : avail_section(11) = "BONDANGLEH"
106 12 : avail_section(12) = "BONDANGLE"
107 12 : avail_section(13) = "IMPDIHEDRALTYPE"
108 12 : avail_section(14) = "IMPDIHEDRALH"
109 12 : avail_section(15) = "IMPDIHEDRAL"
110 12 : avail_section(16) = "DIHEDRALTYPE"
111 12 : avail_section(17) = "DIHEDRALH"
112 12 : avail_section(18) = "DIHEDRAL"
113 12 : avail_section(19) = "LJPARAMETERS"
114 12 : avail_section(20) = "SOLVENTATOM"
115 12 : avail_section(21) = "SOLVENTCONSTR"
116 :
117 12 : gro_info => ff_type%gro_info
118 12 : gro_info%ff_gromos_type = ff_type%ff_type
119 12 : NULLIFY (namearray)
120 : ! ATOMTYPENAME SECTION
121 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the ATOMTYPENAME section'
122 12 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
123 12 : label = TRIM(avail_section(4))
124 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
125 12 : IF (found) THEN
126 12 : CALL parser_get_next_line(parser, 1)
127 12 : CALL parser_get_object(parser, ntype)
128 12 : CALL reallocate(namearray, 1, ntype)
129 96 : DO itype = 1, ntype
130 84 : CALL parser_get_next_line(parser, 1)
131 84 : CALL parser_get_object(parser, namearray(itype), lower_to_upper=.TRUE.)
132 96 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(namearray(itype))
133 : END DO
134 : END IF
135 12 : CALL parser_release(parser)
136 :
137 : ! SOLVENTCONSTR SECTION
138 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the SOLVENTATOM section'
139 12 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
140 :
141 12 : label = TRIM(avail_section(21))
142 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
143 12 : IF (found) THEN
144 8 : CALL parser_get_next_line(parser, 1)
145 8 : CALL parser_get_object(parser, ncon)
146 8 : CALL reallocate(gro_info%solvent_k, 1, ncon)
147 8 : CALL reallocate(gro_info%solvent_r0, 1, ncon)
148 32 : DO icon = 1, ncon
149 24 : CALL parser_get_next_line(parser, 1)
150 24 : CALL parser_get_object(parser, itemp)
151 24 : CALL parser_get_object(parser, itemp)
152 24 : CALL parser_get_object(parser, gro_info%solvent_r0(icon))
153 56 : gro_info%solvent_k(icon) = 0.0_dp
154 : END DO
155 : END IF
156 12 : CALL parser_release(parser)
157 :
158 12 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
159 : ! BONDTYPE SECTION
160 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the BONDTYPE section'
161 12 : label = TRIM(avail_section(7))
162 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
163 12 : IF (found) THEN
164 12 : CALL parser_get_next_line(parser, 1)
165 12 : CALL parser_get_object(parser, ntype)
166 12 : CALL reallocate(gro_info%bond_k, 1, ntype)
167 12 : CALL reallocate(gro_info%bond_r0, 1, ntype)
168 80 : DO itype = 1, ntype
169 68 : CALL parser_get_next_line(parser, 1)
170 68 : CALL parser_get_object(parser, gro_info%bond_k(itype))
171 68 : CALL parser_get_object(parser, gro_info%bond_r0(itype))
172 68 : IF (ff_type%ff_type == do_ff_g96) THEN
173 36 : gro_info%bond_k(itype) = cp_unit_to_cp2k(gro_info%bond_k(itype), "kjmol*nm^-4")
174 : ELSE ! Assume its G87
175 32 : gro_info%bond_k(itype) = (2.0_dp)*gro_info%bond_k(itype)*gro_info%bond_r0(itype)**2
176 32 : gro_info%bond_k(itype) = cp_unit_to_cp2k(gro_info%bond_k(itype), "kjmol*nm^-2")
177 : END IF
178 68 : gro_info%bond_r0(itype) = cp_unit_to_cp2k(gro_info%bond_r0(itype), "nm")
179 80 : IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT BONDTYPE INFO HERE!"
180 : END DO
181 : END IF
182 :
183 : ! BONDANGLETYPE SECTION
184 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the BONDANGLETYPE section'
185 12 : label = TRIM(avail_section(10))
186 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
187 12 : IF (found) THEN
188 12 : CALL parser_get_next_line(parser, 1)
189 12 : CALL parser_get_object(parser, ntype)
190 12 : CALL reallocate(gro_info%bend_k, 1, ntype)
191 12 : CALL reallocate(gro_info%bend_theta0, 1, ntype)
192 104 : DO itype = 1, ntype
193 92 : CALL parser_get_next_line(parser, 1)
194 92 : CALL parser_get_object(parser, gro_info%bend_k(itype))
195 92 : CALL parser_get_object(parser, gro_info%bend_theta0(itype))
196 92 : gro_info%bend_theta0(itype) = cp_unit_to_cp2k(gro_info%bend_theta0(itype), "deg")
197 92 : IF (ff_type%ff_type == do_ff_g96) THEN
198 48 : gro_info%bend_theta0(itype) = COS(gro_info%bend_theta0(itype))
199 : ELSE ! Assume its G87
200 44 : cost2 = COS(gro_info%bend_theta0(itype))*COS(gro_info%bend_theta0(itype))
201 44 : sdet = cost2*cost2 - (2.0_dp*cost2 - 1.0_dp)*(1.0_dp - ekt/gro_info%bend_k(itype))
202 44 : csq = (cost2 - SQRT(sdet))/(2.0_dp*cost2 - 1.0_dp)
203 44 : gro_info%bend_k(itype) = ekt/ACOS(csq)**2
204 : END IF
205 92 : gro_info%bend_k(itype) = cp_unit_to_cp2k(gro_info%bend_k(itype), "kjmol")
206 104 : IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT BONDANGLETYPE INFO HERE!"
207 : END DO
208 : END IF
209 :
210 : ! IMPDIHEDRALTYPE SECTION
211 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the IMPDIHEDRALTYPE section'
212 12 : label = TRIM(avail_section(13))
213 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
214 12 : IF (found) THEN
215 12 : CALL parser_get_next_line(parser, 1)
216 12 : CALL parser_get_object(parser, ntype)
217 12 : CALL reallocate(gro_info%impr_k, 1, ntype)
218 12 : CALL reallocate(gro_info%impr_phi0, 1, ntype)
219 12 : DO itype = 1, ntype
220 0 : CALL parser_get_next_line(parser, 1)
221 0 : CALL parser_get_object(parser, gro_info%impr_k(itype))
222 0 : CALL parser_get_object(parser, gro_info%impr_phi0(itype))
223 0 : gro_info%impr_phi0(itype) = cp_unit_to_cp2k(gro_info%impr_phi0(itype), "deg")
224 0 : gro_info%impr_k(itype) = cp_unit_to_cp2k(gro_info%impr_k(itype), "kjmol*deg^-2")
225 12 : IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT IMPDIHEDRALTYPE INFO HERE!"
226 : END DO
227 : END IF
228 :
229 : ! DIHEDRALTYPE SECTION
230 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the DIHEDRALTYPE section'
231 12 : label = TRIM(avail_section(16))
232 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
233 12 : IF (found) THEN
234 12 : CALL parser_get_next_line(parser, 1)
235 12 : CALL parser_get_object(parser, ntype)
236 12 : CALL reallocate(gro_info%torsion_k, 1, ntype)
237 12 : CALL reallocate(gro_info%torsion_m, 1, ntype)
238 12 : CALL reallocate(gro_info%torsion_phi0, 1, ntype)
239 164 : DO itype = 1, ntype
240 152 : CALL parser_get_next_line(parser, 1)
241 152 : CALL parser_get_object(parser, gro_info%torsion_k(itype))
242 152 : CALL parser_get_object(parser, cosphi0)
243 152 : CALL parser_get_object(parser, gro_info%torsion_m(itype))
244 152 : gro_info%torsion_phi0(itype) = ACOS(cosphi0)
245 152 : gro_info%torsion_k(itype) = cp_unit_to_cp2k(gro_info%torsion_k(itype), "kjmol")
246 316 : IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT DIHEDRALTYPE INFO HERE!"
247 : END DO
248 : END IF
249 :
250 12 : CALL parser_release(parser)
251 :
252 : ! LJPARAMETERS SECTION
253 12 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the LJPARAMETERS section'
254 12 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
255 12 : label = TRIM(avail_section(19))
256 12 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
257 12 : IF (found) THEN
258 12 : CALL parser_get_next_line(parser, 1)
259 12 : CALL parser_get_object(parser, ntype)
260 12 : offset = 0
261 12 : IF (ASSOCIATED(gro_info%nonbond_a)) offset = SIZE(gro_info%nonbond_a)
262 12 : ntype = SIZE(namearray)
263 12 : CALL reallocate(gro_info%nonbond_a, 1, ntype)
264 12 : CALL reallocate(gro_info%nonbond_a_14, 1, ntype)
265 12 : CALL reallocate(gro_info%nonbond_c6, 1, ntype, 1, ntype)
266 12 : CALL reallocate(gro_info%nonbond_c12, 1, ntype, 1, ntype)
267 12 : CALL reallocate(gro_info%nonbond_c6_14, 1, ntype, 1, ntype)
268 12 : CALL reallocate(gro_info%nonbond_c12_14, 1, ntype, 1, ntype)
269 :
270 780 : gro_info%nonbond_c12 = 0._dp
271 780 : gro_info%nonbond_c6 = 0._dp
272 780 : gro_info%nonbond_c12_14 = 0._dp
273 780 : gro_info%nonbond_c6_14 = 0._dp
274 :
275 396 : DO itype = 1, ntype*(ntype + 1)/2
276 384 : CALL parser_get_next_line(parser, 1)
277 384 : CALL parser_get_object(parser, iatom)
278 384 : CALL parser_get_object(parser, jatom)
279 384 : IF (iatom == jatom) THEN
280 84 : gro_info%nonbond_a(iatom) = namearray(iatom)
281 84 : gro_info%nonbond_a_14(iatom) = namearray(iatom)
282 : END IF
283 384 : CALL parser_get_object(parser, gro_info%nonbond_c12(iatom, jatom))
284 384 : CALL parser_get_object(parser, gro_info%nonbond_c6(iatom, jatom))
285 384 : CALL parser_get_object(parser, gro_info%nonbond_c12_14(iatom, jatom))
286 384 : CALL parser_get_object(parser, gro_info%nonbond_c6_14(iatom, jatom))
287 384 : gro_info%nonbond_c6(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c6(iatom, jatom), "kjmol*nm^6")
288 384 : gro_info%nonbond_c12(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c12(iatom, jatom), "kjmol*nm^12")
289 384 : gro_info%nonbond_c6_14(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c6_14(iatom, jatom), "kjmol*nm^6")
290 384 : gro_info%nonbond_c12_14(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c12_14(iatom, jatom), "kjmol*nm^12")
291 :
292 384 : gro_info%nonbond_c6_14(jatom, iatom) = gro_info%nonbond_c6_14(iatom, jatom)
293 384 : gro_info%nonbond_c12_14(jatom, iatom) = gro_info%nonbond_c12_14(iatom, jatom)
294 384 : gro_info%nonbond_c6(jatom, iatom) = gro_info%nonbond_c6(iatom, jatom)
295 384 : gro_info%nonbond_c12(jatom, iatom) = gro_info%nonbond_c12(iatom, jatom)
296 780 : IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT LJPARAMETERS INFO HERE!"
297 : END DO
298 : END IF
299 12 : CALL parser_release(parser)
300 :
301 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
302 12 : "PRINT%FF_INFO")
303 12 : CALL timestop(handle)
304 :
305 12 : DEALLOCATE (namearray)
306 36 : END SUBROUTINE read_force_field_gromos
307 :
308 : ! **************************************************************************************************
309 : !> \brief Reads the charmm force_field
310 : !> \param ff_type ...
311 : !> \param para_env ...
312 : !> \param mm_section ...
313 : !> \author ikuo
314 : ! **************************************************************************************************
315 1756 : SUBROUTINE read_force_field_charmm(ff_type, para_env, mm_section)
316 :
317 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
318 : TYPE(mp_para_env_type), POINTER :: para_env
319 : TYPE(section_vals_type), POINTER :: mm_section
320 :
321 : CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_charmm'
322 :
323 : CHARACTER(LEN=default_string_length) :: label, string, string2, string3, string4
324 : CHARACTER(LEN=default_string_length), DIMENSION(1) :: bond_section
325 : CHARACTER(LEN=default_string_length), DIMENSION(2) :: angl_section, impr_section, &
326 : nbon_section, thet_section
327 : CHARACTER(LEN=default_string_length), &
328 : DIMENSION(20) :: avail_section
329 : INTEGER :: dummy, handle, ilab, iw, nbend, nbond, &
330 : nimpr, nnonbond, nonfo, ntorsion, nub
331 : LOGICAL :: found
332 : TYPE(charmm_info_type), POINTER :: chm_info
333 : TYPE(cp_logger_type), POINTER :: logger
334 : TYPE(cp_parser_type) :: parser
335 :
336 878 : CALL timeset(routineN, handle)
337 878 : NULLIFY (logger)
338 878 : logger => cp_get_default_logger()
339 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
340 878 : extension=".mmLog")
341 :
342 878 : avail_section(1) = "BOND"; bond_section(1) = avail_section(1)
343 878 : avail_section(11) = "BONDS"
344 878 : avail_section(2) = "ANGL"; angl_section(1) = avail_section(2)
345 878 : avail_section(3) = "THETA"; angl_section(2) = avail_section(3)
346 878 : avail_section(12) = "THETAS"
347 878 : avail_section(13) = "ANGLE"
348 878 : avail_section(14) = "ANGLES"
349 878 : avail_section(4) = "DIHEDRAL"; thet_section(1) = avail_section(4)
350 878 : avail_section(15) = "DIHEDRALS"
351 878 : avail_section(5) = "PHI"; thet_section(2) = avail_section(5)
352 878 : avail_section(6) = "IMPROPER"; impr_section(1) = avail_section(6)
353 878 : avail_section(20) = "IMPROPERS"
354 878 : avail_section(7) = "IMPH"; impr_section(2) = avail_section(7)
355 878 : avail_section(16) = "IMPHI"
356 878 : avail_section(8) = "NONBONDED"; nbon_section(1) = avail_section(8)
357 878 : avail_section(9) = "NBOND"; nbon_section(2) = avail_section(9)
358 878 : avail_section(10) = "HBOND"
359 878 : avail_section(17) = "NBFIX"
360 878 : avail_section(18) = "CMAP"
361 878 : avail_section(19) = "END"
362 :
363 878 : chm_info => ff_type%chm_info
364 : !-----------------------------------------------------------------------------
365 : !-----------------------------------------------------------------------------
366 : ! 1. Read in all the Bonds info from the param file here
367 : ! Vbond = Kb(b-b0)^2
368 : ! UNITS for Kb: [(kcal/mol)/(A^2)] to [Eh/(AU^2)]
369 : !-----------------------------------------------------------------------------
370 878 : nbond = 0
371 1756 : DO ilab = 1, SIZE(bond_section)
372 878 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
373 878 : label = TRIM(bond_section(ilab))
374 : DO
375 3244 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
376 3244 : IF (found) THEN
377 2366 : CALL parser_get_object(parser, string)
378 2366 : IF (INDEX(string, TRIM(label)) /= 1) CYCLE
379 878 : CALL charmm_get_next_line(parser, 1)
380 19876 : DO
381 20754 : CALL parser_get_object(parser, string)
382 20754 : CALL uppercase(string)
383 429688 : IF (ANY(string == avail_section)) EXIT
384 19876 : CALL parser_get_object(parser, string2)
385 19876 : CALL uppercase(string2)
386 19876 : nbond = nbond + 1
387 19876 : CALL reallocate(chm_info%bond_a, 1, nbond)
388 19876 : CALL reallocate(chm_info%bond_b, 1, nbond)
389 19876 : CALL reallocate(chm_info%bond_k, 1, nbond)
390 19876 : CALL reallocate(chm_info%bond_r0, 1, nbond)
391 19876 : chm_info%bond_a(nbond) = string
392 19876 : chm_info%bond_b(nbond) = string2
393 19876 : CALL parser_get_object(parser, chm_info%bond_k(nbond))
394 19876 : CALL parser_get_object(parser, chm_info%bond_r0(nbond))
395 20907 : IF (iw > 0) WRITE (iw, *) " CHM BOND ", nbond, &
396 1031 : TRIM(chm_info%bond_a(nbond)), " ", &
397 1031 : TRIM(chm_info%bond_b(nbond)), " ", &
398 1031 : chm_info%bond_k(nbond), &
399 2062 : chm_info%bond_r0(nbond)
400 : ! Do some units conversion into internal atomic units
401 19876 : chm_info%bond_r0(nbond) = cp_unit_to_cp2k(chm_info%bond_r0(nbond), "angstrom")
402 19876 : chm_info%bond_k(nbond) = cp_unit_to_cp2k(chm_info%bond_k(nbond), "kcalmol*angstrom^-2")
403 22242 : CALL charmm_get_next_line(parser, 1)
404 : END DO
405 : ELSE
406 : EXIT
407 : END IF
408 : END DO
409 2634 : CALL parser_release(parser)
410 : END DO
411 : !-----------------------------------------------------------------------------
412 : !-----------------------------------------------------------------------------
413 : ! 2. Read in all the Bends and UB info from the param file here
414 : ! Vangle = Ktheta(theta-theta0)^2
415 : ! UNITS for Ktheta: [(kcal/mol)/(rad^2)] to [Eh/(rad^2)]
416 : ! FACTOR of "2" rolled into Ktheta
417 : ! Vub = Kub(S-S0)^2
418 : ! UNITS for Kub: [(kcal/mol)/(A^2)] to [Eh/(AU^2)]
419 : !-----------------------------------------------------------------------------
420 878 : nbend = 0
421 878 : nub = 0
422 2634 : DO ilab = 1, SIZE(angl_section)
423 1756 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
424 1756 : label = TRIM(angl_section(ilab))
425 : DO
426 2644 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
427 2644 : IF (found) THEN
428 888 : CALL parser_get_object(parser, string)
429 888 : IF (INDEX(string, TRIM(label)) /= 1) CYCLE
430 878 : CALL charmm_get_next_line(parser, 1)
431 : DO
432 45533 : CALL parser_get_object(parser, string)
433 45533 : CALL uppercase(string)
434 950825 : IF (ANY(string == avail_section)) EXIT
435 44655 : CALL parser_get_object(parser, string2)
436 44655 : CALL parser_get_object(parser, string3)
437 44655 : CALL uppercase(string2)
438 44655 : CALL uppercase(string3)
439 44655 : nbend = nbend + 1
440 44655 : CALL reallocate(chm_info%bend_a, 1, nbend)
441 44655 : CALL reallocate(chm_info%bend_b, 1, nbend)
442 44655 : CALL reallocate(chm_info%bend_c, 1, nbend)
443 44655 : CALL reallocate(chm_info%bend_k, 1, nbend)
444 44655 : CALL reallocate(chm_info%bend_theta0, 1, nbend)
445 44655 : chm_info%bend_a(nbend) = string
446 44655 : chm_info%bend_b(nbend) = string2
447 44655 : chm_info%bend_c(nbend) = string3
448 44655 : CALL parser_get_object(parser, chm_info%bend_k(nbend))
449 44655 : CALL parser_get_object(parser, chm_info%bend_theta0(nbend))
450 46630 : IF (iw > 0) WRITE (iw, *) " CHM BEND ", nbend, &
451 1975 : TRIM(chm_info%bend_a(nbend)), " ", &
452 1975 : TRIM(chm_info%bend_b(nbend)), " ", &
453 1975 : TRIM(chm_info%bend_c(nbend)), " ", &
454 1975 : chm_info%bend_k(nbend), &
455 3950 : chm_info%bend_theta0(nbend)
456 : ! Do some units conversion into internal atomic units
457 44655 : chm_info%bend_theta0(nbend) = cp_unit_to_cp2k(chm_info%bend_theta0(nbend), "deg")
458 44655 : chm_info%bend_k(nbend) = cp_unit_to_cp2k(chm_info%bend_k(nbend), "kcalmol*rad^-2")
459 44655 : IF (parser_test_next_token(parser) == "FLT") THEN
460 13074 : nub = nub + 1
461 13074 : CALL reallocate(chm_info%ub_a, 1, nub)
462 13074 : CALL reallocate(chm_info%ub_b, 1, nub)
463 13074 : CALL reallocate(chm_info%ub_c, 1, nub)
464 13074 : CALL reallocate(chm_info%ub_k, 1, nub)
465 13074 : CALL reallocate(chm_info%ub_r0, 1, nub)
466 13074 : chm_info%ub_a(nub) = string
467 13074 : chm_info%ub_b(nub) = string2
468 13074 : chm_info%ub_c(nub) = string3
469 13074 : CALL parser_get_object(parser, chm_info%ub_k(nub))
470 13074 : CALL parser_get_object(parser, chm_info%ub_r0(nub))
471 13390 : IF (iw > 0) WRITE (iw, *) " CHM UB ", nub, &
472 316 : TRIM(chm_info%ub_a(nub)), " ", &
473 316 : TRIM(chm_info%ub_b(nub)), " ", &
474 316 : TRIM(chm_info%ub_c(nub)), " ", &
475 316 : chm_info%ub_k(nub), &
476 632 : chm_info%ub_r0(nub)
477 : ! Do some units conversion into internal atomic units
478 13074 : chm_info%ub_r0(nub) = cp_unit_to_cp2k(chm_info%ub_r0(nub), "angstrom")
479 57729 : chm_info%ub_k(nub) = cp_unit_to_cp2k(chm_info%ub_k(nub), "kcalmol*angstrom^-2")
480 : END IF
481 45543 : CALL charmm_get_next_line(parser, 1)
482 : END DO
483 : ELSE
484 : EXIT
485 : END IF
486 : END DO
487 4390 : CALL parser_release(parser)
488 : END DO
489 : !-----------------------------------------------------------------------------
490 : !-----------------------------------------------------------------------------
491 : ! 3. Read in all the Dihedrals info from the param file here
492 : ! Vtorsion = Kphi(1+COS(n(phi)-delta))
493 : ! UNITS for Kphi: [(kcal/mol)] to [Eh]
494 : !-----------------------------------------------------------------------------
495 878 : ntorsion = 0
496 2634 : DO ilab = 1, SIZE(thet_section)
497 1756 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
498 1756 : label = TRIM(thet_section(ilab))
499 : DO
500 2644 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
501 2644 : IF (found) THEN
502 888 : CALL parser_get_object(parser, string)
503 888 : IF (INDEX(string, TRIM(label)) /= 1) CYCLE
504 878 : CALL charmm_get_next_line(parser, 1)
505 52528 : DO
506 53406 : CALL parser_get_object(parser, string)
507 53406 : CALL uppercase(string)
508 1108456 : IF (ANY(string == avail_section)) EXIT
509 52528 : CALL parser_get_object(parser, string2)
510 52528 : CALL parser_get_object(parser, string3)
511 52528 : CALL parser_get_object(parser, string4)
512 52528 : CALL uppercase(string2)
513 52528 : CALL uppercase(string3)
514 52528 : CALL uppercase(string4)
515 52528 : ntorsion = ntorsion + 1
516 52528 : CALL reallocate(chm_info%torsion_a, 1, ntorsion)
517 52528 : CALL reallocate(chm_info%torsion_b, 1, ntorsion)
518 52528 : CALL reallocate(chm_info%torsion_c, 1, ntorsion)
519 52528 : CALL reallocate(chm_info%torsion_d, 1, ntorsion)
520 52528 : CALL reallocate(chm_info%torsion_k, 1, ntorsion)
521 52528 : CALL reallocate(chm_info%torsion_m, 1, ntorsion)
522 52528 : CALL reallocate(chm_info%torsion_phi0, 1, ntorsion)
523 52528 : chm_info%torsion_a(ntorsion) = string
524 52528 : chm_info%torsion_b(ntorsion) = string2
525 52528 : chm_info%torsion_c(ntorsion) = string3
526 52528 : chm_info%torsion_d(ntorsion) = string4
527 52528 : CALL parser_get_object(parser, chm_info%torsion_k(ntorsion))
528 52528 : CALL parser_get_object(parser, chm_info%torsion_m(ntorsion))
529 52528 : CALL parser_get_object(parser, chm_info%torsion_phi0(ntorsion))
530 54922 : IF (iw > 0) WRITE (iw, *) " CHM TORSION ", ntorsion, &
531 2394 : TRIM(chm_info%torsion_a(ntorsion)), " ", &
532 2394 : TRIM(chm_info%torsion_b(ntorsion)), " ", &
533 2394 : TRIM(chm_info%torsion_c(ntorsion)), " ", &
534 2394 : TRIM(chm_info%torsion_d(ntorsion)), " ", &
535 2394 : chm_info%torsion_k(ntorsion), &
536 2394 : chm_info%torsion_m(ntorsion), &
537 4788 : chm_info%torsion_phi0(ntorsion)
538 : ! Do some units conversion into internal atomic units
539 : chm_info%torsion_phi0(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_phi0(ntorsion), &
540 52528 : "deg")
541 52528 : chm_info%torsion_k(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_k(ntorsion), "kcalmol")
542 53416 : CALL charmm_get_next_line(parser, 1)
543 : END DO
544 : ELSE
545 : EXIT
546 : END IF
547 : END DO
548 4390 : CALL parser_release(parser)
549 : END DO
550 : !-----------------------------------------------------------------------------
551 : !-----------------------------------------------------------------------------
552 : ! 4. Read in all the Improper info from the param file here
553 : ! Vimpr = Kpsi(psi-psi0)^2
554 : ! UNITS for Kpsi: [(kcal/mol)/(rad^2)] to [Eh/(rad^2)]
555 : !-----------------------------------------------------------------------------
556 878 : nimpr = 0
557 2634 : DO ilab = 1, SIZE(impr_section)
558 1756 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
559 1756 : label = TRIM(impr_section(ilab))
560 : DO
561 2634 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
562 2634 : IF (found) THEN
563 878 : CALL parser_get_object(parser, string)
564 878 : IF (INDEX(string, TRIM(label)) /= 1) CYCLE
565 878 : CALL charmm_get_next_line(parser, 1)
566 4156 : DO
567 5034 : CALL parser_get_object(parser, string)
568 5034 : CALL uppercase(string)
569 94300 : IF (ANY(string == avail_section)) EXIT
570 4156 : CALL parser_get_object(parser, string2)
571 4156 : CALL parser_get_object(parser, string3)
572 4156 : CALL parser_get_object(parser, string4)
573 4156 : CALL uppercase(string2)
574 4156 : CALL uppercase(string3)
575 4156 : CALL uppercase(string4)
576 4156 : nimpr = nimpr + 1
577 4156 : CALL reallocate(chm_info%impr_a, 1, nimpr)
578 4156 : CALL reallocate(chm_info%impr_b, 1, nimpr)
579 4156 : CALL reallocate(chm_info%impr_c, 1, nimpr)
580 4156 : CALL reallocate(chm_info%impr_d, 1, nimpr)
581 4156 : CALL reallocate(chm_info%impr_k, 1, nimpr)
582 4156 : CALL reallocate(chm_info%impr_phi0, 1, nimpr)
583 4156 : chm_info%impr_a(nimpr) = string
584 4156 : chm_info%impr_b(nimpr) = string2
585 4156 : chm_info%impr_c(nimpr) = string3
586 4156 : chm_info%impr_d(nimpr) = string4
587 4156 : CALL parser_get_object(parser, chm_info%impr_k(nimpr))
588 4156 : CALL parser_get_object(parser, dummy)
589 4156 : CALL parser_get_object(parser, chm_info%impr_phi0(nimpr))
590 4262 : IF (iw > 0) WRITE (iw, *) " CHM IMPROPERS ", nimpr, &
591 106 : TRIM(chm_info%impr_a(nimpr)), " ", &
592 106 : TRIM(chm_info%impr_b(nimpr)), " ", &
593 106 : TRIM(chm_info%impr_c(nimpr)), " ", &
594 106 : TRIM(chm_info%impr_d(nimpr)), " ", &
595 106 : chm_info%impr_k(nimpr), &
596 212 : chm_info%impr_phi0(nimpr)
597 : ! Do some units conversion into internal atomic units
598 4156 : chm_info%impr_phi0(nimpr) = cp_unit_to_cp2k(chm_info%impr_phi0(nimpr), "deg")
599 4156 : chm_info%impr_k(nimpr) = cp_unit_to_cp2k(chm_info%impr_k(nimpr), "kcalmol")
600 5034 : CALL charmm_get_next_line(parser, 1)
601 : END DO
602 : ELSE
603 : EXIT
604 : END IF
605 : END DO
606 4390 : CALL parser_release(parser)
607 : END DO
608 : !-----------------------------------------------------------------------------
609 : !-----------------------------------------------------------------------------
610 : ! 5. Read in all the Nonbonded info from the param file here
611 : !-----------------------------------------------------------------------------
612 878 : nnonbond = 0
613 878 : nonfo = 0
614 2634 : DO ilab = 1, SIZE(nbon_section)
615 1756 : CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
616 1756 : label = TRIM(nbon_section(ilab))
617 : DO
618 3512 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
619 3512 : IF (found) THEN
620 1756 : CALL parser_get_object(parser, string)
621 1756 : IF (INDEX(string, TRIM(label)) /= 1) CYCLE
622 878 : CALL charmm_get_next_line(parser, 1)
623 : DO
624 12714 : CALL parser_get_object(parser, string)
625 12714 : CALL uppercase(string)
626 259928 : IF (ANY(string == avail_section)) EXIT
627 11836 : nnonbond = nnonbond + 1
628 11836 : CALL reallocate(chm_info%nonbond_a, 1, nnonbond)
629 11836 : CALL reallocate(chm_info%nonbond_eps, 1, nnonbond)
630 11836 : CALL reallocate(chm_info%nonbond_rmin2, 1, nnonbond)
631 11836 : chm_info%nonbond_a(nnonbond) = string
632 11836 : CALL parser_get_object(parser, chm_info%nonbond_eps(nnonbond))
633 11836 : CALL parser_get_object(parser, chm_info%nonbond_eps(nnonbond))
634 11836 : CALL parser_get_object(parser, chm_info%nonbond_rmin2(nnonbond))
635 12737 : IF (iw > 0) WRITE (iw, *) " CHM NONBOND ", nnonbond, &
636 901 : TRIM(chm_info%nonbond_a(nnonbond)), " ", &
637 901 : chm_info%nonbond_eps(nnonbond), &
638 1802 : chm_info%nonbond_rmin2(nnonbond)
639 : chm_info%nonbond_rmin2(nnonbond) = cp_unit_to_cp2k(chm_info%nonbond_rmin2(nnonbond), &
640 11836 : "angstrom")
641 : chm_info%nonbond_eps(nnonbond) = cp_unit_to_cp2k(chm_info%nonbond_eps(nnonbond), &
642 11836 : "kcalmol")
643 11836 : IF (parser_test_next_token(parser) == "FLT") THEN
644 2198 : nonfo = nonfo + 1
645 2198 : CALL reallocate(chm_info%nonbond_a_14, 1, nonfo)
646 2198 : CALL reallocate(chm_info%nonbond_eps_14, 1, nonfo)
647 2198 : CALL reallocate(chm_info%nonbond_rmin2_14, 1, nonfo)
648 2198 : chm_info%nonbond_a_14(nonfo) = chm_info%nonbond_a(nnonbond)
649 2198 : CALL parser_get_object(parser, chm_info%nonbond_eps_14(nonfo))
650 2198 : CALL parser_get_object(parser, chm_info%nonbond_eps_14(nonfo))
651 2198 : CALL parser_get_object(parser, chm_info%nonbond_rmin2_14(nonfo))
652 2242 : IF (iw > 0) WRITE (iw, *) " CHM ONFO ", nonfo, &
653 44 : TRIM(chm_info%nonbond_a_14(nonfo)), " ", &
654 44 : chm_info%nonbond_eps_14(nonfo), &
655 88 : chm_info%nonbond_rmin2_14(nonfo)
656 : chm_info%nonbond_rmin2_14(nonfo) = cp_unit_to_cp2k(chm_info%nonbond_rmin2_14(nonfo), &
657 2198 : "angstrom")
658 : chm_info%nonbond_eps_14(nonfo) = cp_unit_to_cp2k(chm_info%nonbond_eps_14(nonfo), &
659 14034 : "kcalmol")
660 : END IF
661 13592 : CALL charmm_get_next_line(parser, 1)
662 : END DO
663 : ELSE
664 : EXIT
665 : END IF
666 : END DO
667 4390 : CALL parser_release(parser)
668 : END DO
669 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
670 878 : "PRINT%FF_INFO")
671 878 : CALL timestop(handle)
672 :
673 2634 : END SUBROUTINE read_force_field_charmm
674 :
675 : ! **************************************************************************************************
676 : !> \brief Reads the AMBER force_field
677 : !> \param ff_type ...
678 : !> \param para_env ...
679 : !> \param mm_section ...
680 : !> \param particle_set ...
681 : !> \author Teodoro Laino [tlaino, teodoro.laino-AT-gmail.com] - 11.2008
682 : ! **************************************************************************************************
683 14 : SUBROUTINE read_force_field_amber(ff_type, para_env, mm_section, particle_set)
684 :
685 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
686 : TYPE(mp_para_env_type), POINTER :: para_env
687 : TYPE(section_vals_type), POINTER :: mm_section
688 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
689 :
690 : CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_amber'
691 :
692 : INTEGER :: handle, i, iw
693 : TYPE(amber_info_type), POINTER :: amb_info
694 : TYPE(cp_logger_type), POINTER :: logger
695 :
696 14 : CALL timeset(routineN, handle)
697 14 : NULLIFY (logger)
698 14 : logger => cp_get_default_logger()
699 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
700 14 : extension=".mmLog")
701 :
702 14 : amb_info => ff_type%amb_info
703 :
704 : ! Read the Amber topology file to gather the forcefield information
705 : CALL rdparm_amber_8(ff_type%ff_file_name, iw, para_env, do_connectivity=.FALSE., &
706 14 : do_forcefield=.TRUE., particle_set=particle_set, amb_info=amb_info)
707 :
708 : !-----------------------------------------------------------------------------
709 : ! 1. Converts all the Bonds info from the param file here
710 : ! Vbond = Kb(b-b0)^2
711 : ! UNITS for Kb: [(kcal/mol)/(A^2)] to [Eh/(AU^2)]
712 : !-----------------------------------------------------------------------------
713 1738 : DO i = 1, SIZE(amb_info%bond_a)
714 1737 : IF (iw > 0) WRITE (iw, *) " AMB BOND ", i, &
715 13 : TRIM(amb_info%bond_a(i)), " ", &
716 13 : TRIM(amb_info%bond_b(i)), " ", &
717 13 : amb_info%bond_k(i), &
718 26 : amb_info%bond_r0(i)
719 :
720 : ! Do some units conversion into internal atomic units
721 1724 : amb_info%bond_r0(i) = cp_unit_to_cp2k(amb_info%bond_r0(i), "angstrom")
722 1738 : amb_info%bond_k(i) = cp_unit_to_cp2k(amb_info%bond_k(i), "kcalmol*angstrom^-2")
723 : END DO
724 :
725 : !-----------------------------------------------------------------------------
726 : ! 2. Converts all the Bends info from the param file here
727 : ! Vangle = Ktheta(theta-theta0)^2
728 : ! UNITS for Ktheta: [(kcal/mol)/(rad^2)] to [Eh/(rad^2)]
729 : ! FACTOR of "2" rolled into Ktheta
730 : ! Vub = Kub(S-S0)^2
731 : ! UNITS for Kub: [(kcal/mol)/(A^2)] to [Eh/(AU^2)]
732 : !-----------------------------------------------------------------------------
733 3650 : DO i = 1, SIZE(amb_info%bend_a)
734 3658 : IF (iw > 0) WRITE (iw, *) " AMB BEND ", i, &
735 22 : TRIM(amb_info%bend_a(i)), " ", &
736 22 : TRIM(amb_info%bend_b(i)), " ", &
737 22 : TRIM(amb_info%bend_c(i)), " ", &
738 22 : amb_info%bend_k(i), &
739 44 : amb_info%bend_theta0(i)
740 :
741 : ! Do some units conversion into internal atomic units
742 3636 : amb_info%bend_theta0(i) = cp_unit_to_cp2k(amb_info%bend_theta0(i), "rad")
743 3650 : amb_info%bend_k(i) = cp_unit_to_cp2k(amb_info%bend_k(i), "kcalmol*rad^-2")
744 : END DO
745 :
746 : !-----------------------------------------------------------------------------
747 : ! 3. Converts all the Dihedrals info from the param file here
748 : ! Vtorsion = Kphi(1+COS(n(phi)-delta))
749 : ! UNITS for Kphi: [(kcal/mol)] to [Eh]
750 : !-----------------------------------------------------------------------------
751 10096 : DO i = 1, SIZE(amb_info%torsion_a)
752 10121 : IF (iw > 0) WRITE (iw, *) " AMB TORSION ", i, &
753 39 : TRIM(amb_info%torsion_a(i)), " ", &
754 39 : TRIM(amb_info%torsion_b(i)), " ", &
755 39 : TRIM(amb_info%torsion_c(i)), " ", &
756 39 : TRIM(amb_info%torsion_d(i)), " ", &
757 39 : amb_info%torsion_k(i), &
758 39 : amb_info%torsion_m(i), &
759 78 : amb_info%torsion_phi0(i)
760 :
761 : ! Do some units conversion into internal atomic units
762 10082 : amb_info%torsion_phi0(i) = cp_unit_to_cp2k(amb_info%torsion_phi0(i), "rad")
763 10096 : amb_info%torsion_k(i) = cp_unit_to_cp2k(amb_info%torsion_k(i), "kcalmol")
764 : END DO
765 :
766 : ! Do some units conversion into internal atomic units
767 14 : IF (ASSOCIATED(amb_info%raw_torsion_phi0)) THEN
768 334 : DO i = 1, SIZE(amb_info%raw_torsion_k)
769 322 : amb_info%raw_torsion_phi0(i) = cp_unit_to_cp2k(amb_info%raw_torsion_phi0(i), "rad")
770 334 : amb_info%raw_torsion_k(i) = cp_unit_to_cp2k(amb_info%raw_torsion_k(i), "kcalmol")
771 : END DO
772 : END IF
773 :
774 : !-----------------------------------------------------------------------------
775 : ! 4. Converts all the Nonbonded info from the param file here
776 : !-----------------------------------------------------------------------------
777 1474 : DO i = 1, SIZE(amb_info%nonbond_eps)
778 1472 : IF (iw > 0) WRITE (iw, *) " AMB NONBOND ", i, &
779 12 : TRIM(amb_info%nonbond_a(i)), " ", &
780 12 : amb_info%nonbond_eps(i), &
781 24 : amb_info%nonbond_rmin2(i)
782 :
783 : ! Do some units conversion into internal atomic units
784 1460 : amb_info%nonbond_rmin2(i) = cp_unit_to_cp2k(amb_info%nonbond_rmin2(i), "angstrom")
785 1474 : amb_info%nonbond_eps(i) = cp_unit_to_cp2k(amb_info%nonbond_eps(i), "kcalmol")
786 : END DO
787 14 : CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
788 14 : CALL timestop(handle)
789 14 : END SUBROUTINE read_force_field_amber
790 :
791 : ! **************************************************************************************************
792 : !> \brief This function is simply a wrap to the parser_get_next_line..
793 : !> Comments: This routine would not be necessary if the continuation
794 : !> char for CHARMM would not be the "-".. How can you choose this
795 : !> character in a file of numbers as a continuation char????
796 : !> This sounds simply crazy....
797 : !> \param parser ...
798 : !> \param nline ...
799 : !> \author Teodoro Laino - Zurich University - 06.2007
800 : ! **************************************************************************************************
801 137441 : SUBROUTINE charmm_get_next_line(parser, nline)
802 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
803 : INTEGER, INTENT(IN) :: nline
804 :
805 : CHARACTER(LEN=1), PARAMETER :: continuation_char = "-"
806 :
807 : INTEGER :: i, len_line
808 :
809 274882 : DO i = 1, nline
810 137441 : len_line = LEN_TRIM(parser%input_line)
811 137451 : DO WHILE (parser%input_line(len_line:len_line) == continuation_char)
812 10 : CALL parser_get_next_line(parser, 1)
813 10 : len_line = LEN_TRIM(parser%input_line)
814 : END DO
815 274882 : CALL parser_get_next_line(parser, 1)
816 : END DO
817 :
818 137441 : END SUBROUTINE charmm_get_next_line
819 :
820 : END MODULE force_fields_ext
|
