(git:374b731)
Loading...
Searching...
No Matches
force_fields_ext.F
Go to the documentation of this file.
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! **************************************************************************************************
34 USE cp_units, ONLY: cp_unit_to_cp2k
41 USE kinds, ONLY: default_string_length,&
42 dp
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, &
58
59CONTAINS
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 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 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 CALL timeset(routinen, handle)
90 NULLIFY (logger)
91 logger => cp_get_default_logger()
92 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
93 extension=".mmLog")
94
95 avail_section(1) = "TITLE"
96 avail_section(2) = "TOPPHYSCON"
97 avail_section(3) = "TOPVERSION"
98 avail_section(4) = "ATOMTYPENAME"
99 avail_section(5) = "RESNAME"
100 avail_section(6) = "SOLUTEATOM"
101 avail_section(7) = "BONDTYPE"
102 avail_section(8) = "BONDH"
103 avail_section(9) = "BOND"
104 avail_section(10) = "BONDANGLETYPE"
105 avail_section(11) = "BONDANGLEH"
106 avail_section(12) = "BONDANGLE"
107 avail_section(13) = "IMPDIHEDRALTYPE"
108 avail_section(14) = "IMPDIHEDRALH"
109 avail_section(15) = "IMPDIHEDRAL"
110 avail_section(16) = "DIHEDRALTYPE"
111 avail_section(17) = "DIHEDRALH"
112 avail_section(18) = "DIHEDRAL"
113 avail_section(19) = "LJPARAMETERS"
114 avail_section(20) = "SOLVENTATOM"
115 avail_section(21) = "SOLVENTCONSTR"
116
117 gro_info => ff_type%gro_info
118 gro_info%ff_gromos_type = ff_type%ff_type
119 NULLIFY (namearray)
120 ! ATOMTYPENAME SECTION
121 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the ATOMTYPENAME section'
122 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
123 label = trim(avail_section(4))
124 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
125 IF (found) THEN
126 CALL parser_get_next_line(parser, 1)
127 CALL parser_get_object(parser, ntype)
128 CALL reallocate(namearray, 1, ntype)
129 DO itype = 1, ntype
130 CALL parser_get_next_line(parser, 1)
131 CALL parser_get_object(parser, namearray(itype), lower_to_upper=.true.)
132 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", trim(namearray(itype))
133 END DO
134 END IF
135 CALL parser_release(parser)
136
137 ! SOLVENTCONSTR SECTION
138 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the SOLVENTATOM section'
139 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
140
141 label = trim(avail_section(21))
142 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
143 IF (found) THEN
144 CALL parser_get_next_line(parser, 1)
145 CALL parser_get_object(parser, ncon)
146 CALL reallocate(gro_info%solvent_k, 1, ncon)
147 CALL reallocate(gro_info%solvent_r0, 1, ncon)
148 DO icon = 1, ncon
149 CALL parser_get_next_line(parser, 1)
150 CALL parser_get_object(parser, itemp)
151 CALL parser_get_object(parser, itemp)
152 CALL parser_get_object(parser, gro_info%solvent_r0(icon))
153 gro_info%solvent_k(icon) = 0.0_dp
154 END DO
155 END IF
156 CALL parser_release(parser)
157
158 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
159 ! BONDTYPE SECTION
160 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the BONDTYPE section'
161 label = trim(avail_section(7))
162 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
163 IF (found) THEN
164 CALL parser_get_next_line(parser, 1)
165 CALL parser_get_object(parser, ntype)
166 CALL reallocate(gro_info%bond_k, 1, ntype)
167 CALL reallocate(gro_info%bond_r0, 1, ntype)
168 DO itype = 1, ntype
169 CALL parser_get_next_line(parser, 1)
170 CALL parser_get_object(parser, gro_info%bond_k(itype))
171 CALL parser_get_object(parser, gro_info%bond_r0(itype))
172 IF (ff_type%ff_type == do_ff_g96) THEN
173 gro_info%bond_k(itype) = cp_unit_to_cp2k(gro_info%bond_k(itype), "kjmol*nm^-4")
174 ELSE ! Assume its G87
175 gro_info%bond_k(itype) = (2.0_dp)*gro_info%bond_k(itype)*gro_info%bond_r0(itype)**2
176 gro_info%bond_k(itype) = cp_unit_to_cp2k(gro_info%bond_k(itype), "kjmol*nm^-2")
177 END IF
178 gro_info%bond_r0(itype) = cp_unit_to_cp2k(gro_info%bond_r0(itype), "nm")
179 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT BONDTYPE INFO HERE!!!!"
180 END DO
181 END IF
182
183 ! BONDANGLETYPE SECTION
184 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the BONDANGLETYPE section'
185 label = trim(avail_section(10))
186 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
187 IF (found) THEN
188 CALL parser_get_next_line(parser, 1)
189 CALL parser_get_object(parser, ntype)
190 CALL reallocate(gro_info%bend_k, 1, ntype)
191 CALL reallocate(gro_info%bend_theta0, 1, ntype)
192 DO itype = 1, ntype
193 CALL parser_get_next_line(parser, 1)
194 CALL parser_get_object(parser, gro_info%bend_k(itype))
195 CALL parser_get_object(parser, gro_info%bend_theta0(itype))
196 gro_info%bend_theta0(itype) = cp_unit_to_cp2k(gro_info%bend_theta0(itype), "deg")
197 IF (ff_type%ff_type == do_ff_g96) THEN
198 gro_info%bend_theta0(itype) = cos(gro_info%bend_theta0(itype))
199 ELSE ! Assume its G87
200 cost2 = cos(gro_info%bend_theta0(itype))*cos(gro_info%bend_theta0(itype))
201 sdet = cost2*cost2 - (2.0_dp*cost2 - 1.0_dp)*(1.0_dp - ekt/gro_info%bend_k(itype))
202 csq = (cost2 - sqrt(sdet))/(2.0_dp*cost2 - 1.0_dp)
203 gro_info%bend_k(itype) = ekt/acos(csq)**2
204 END IF
205 gro_info%bend_k(itype) = cp_unit_to_cp2k(gro_info%bend_k(itype), "kjmol")
206 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT BONDANGLETYPE INFO HERE!!!!"
207 END DO
208 END IF
209
210 ! IMPDIHEDRALTYPE SECTION
211 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the IMPDIHEDRALTYPE section'
212 label = trim(avail_section(13))
213 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
214 IF (found) THEN
215 CALL parser_get_next_line(parser, 1)
216 CALL parser_get_object(parser, ntype)
217 CALL reallocate(gro_info%impr_k, 1, ntype)
218 CALL reallocate(gro_info%impr_phi0, 1, ntype)
219 DO itype = 1, ntype
220 CALL parser_get_next_line(parser, 1)
221 CALL parser_get_object(parser, gro_info%impr_k(itype))
222 CALL parser_get_object(parser, gro_info%impr_phi0(itype))
223 gro_info%impr_phi0(itype) = cp_unit_to_cp2k(gro_info%impr_phi0(itype), "deg")
224 gro_info%impr_k(itype) = cp_unit_to_cp2k(gro_info%impr_k(itype), "kjmol*deg^-2")
225 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT IMPDIHEDRALTYPE INFO HERE!!!!"
226 END DO
227 END IF
228
229 ! DIHEDRALTYPE SECTION
230 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the DIHEDRALTYPE section'
231 label = trim(avail_section(16))
232 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
233 IF (found) THEN
234 CALL parser_get_next_line(parser, 1)
235 CALL parser_get_object(parser, ntype)
236 CALL reallocate(gro_info%torsion_k, 1, ntype)
237 CALL reallocate(gro_info%torsion_m, 1, ntype)
238 CALL reallocate(gro_info%torsion_phi0, 1, ntype)
239 DO itype = 1, ntype
240 CALL parser_get_next_line(parser, 1)
241 CALL parser_get_object(parser, gro_info%torsion_k(itype))
242 CALL parser_get_object(parser, cosphi0)
243 CALL parser_get_object(parser, gro_info%torsion_m(itype))
244 gro_info%torsion_phi0(itype) = acos(cosphi0)
245 gro_info%torsion_k(itype) = cp_unit_to_cp2k(gro_info%torsion_k(itype), "kjmol")
246 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT DIHEDRALTYPE INFO HERE!!!!"
247 END DO
248 END IF
249
250 CALL parser_release(parser)
251
252 ! LJPARAMETERS SECTION
253 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the LJPARAMETERS section'
254 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
255 label = trim(avail_section(19))
256 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
257 IF (found) THEN
258 CALL parser_get_next_line(parser, 1)
259 CALL parser_get_object(parser, ntype)
260 offset = 0
261 IF (ASSOCIATED(gro_info%nonbond_a)) offset = SIZE(gro_info%nonbond_a)
262 ntype = SIZE(namearray)
263 CALL reallocate(gro_info%nonbond_a, 1, ntype)
264 CALL reallocate(gro_info%nonbond_a_14, 1, ntype)
265 CALL reallocate(gro_info%nonbond_c6, 1, ntype, 1, ntype)
266 CALL reallocate(gro_info%nonbond_c12, 1, ntype, 1, ntype)
267 CALL reallocate(gro_info%nonbond_c6_14, 1, ntype, 1, ntype)
268 CALL reallocate(gro_info%nonbond_c12_14, 1, ntype, 1, ntype)
269
270 gro_info%nonbond_c12 = 0._dp
271 gro_info%nonbond_c6 = 0._dp
272 gro_info%nonbond_c12_14 = 0._dp
273 gro_info%nonbond_c6_14 = 0._dp
274
275 DO itype = 1, ntype*(ntype + 1)/2
276 CALL parser_get_next_line(parser, 1)
277 CALL parser_get_object(parser, iatom)
278 CALL parser_get_object(parser, jatom)
279 IF (iatom == jatom) THEN
280 gro_info%nonbond_a(iatom) = namearray(iatom)
281 gro_info%nonbond_a_14(iatom) = namearray(iatom)
282 END IF
283 CALL parser_get_object(parser, gro_info%nonbond_c12(iatom, jatom))
284 CALL parser_get_object(parser, gro_info%nonbond_c6(iatom, jatom))
285 CALL parser_get_object(parser, gro_info%nonbond_c12_14(iatom, jatom))
286 CALL parser_get_object(parser, gro_info%nonbond_c6_14(iatom, jatom))
287 gro_info%nonbond_c6(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c6(iatom, jatom), "kjmol*nm^6")
288 gro_info%nonbond_c12(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c12(iatom, jatom), "kjmol*nm^12")
289 gro_info%nonbond_c6_14(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c6_14(iatom, jatom), "kjmol*nm^6")
290 gro_info%nonbond_c12_14(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c12_14(iatom, jatom), "kjmol*nm^12")
291
292 gro_info%nonbond_c6_14(jatom, iatom) = gro_info%nonbond_c6_14(iatom, jatom)
293 gro_info%nonbond_c12_14(jatom, iatom) = gro_info%nonbond_c12_14(iatom, jatom)
294 gro_info%nonbond_c6(jatom, iatom) = gro_info%nonbond_c6(iatom, jatom)
295 gro_info%nonbond_c12(jatom, iatom) = gro_info%nonbond_c12(iatom, jatom)
296 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT LJPARAMETERS INFO HERE!!!!"
297 END DO
298 END IF
299 CALL parser_release(parser)
300
301 CALL cp_print_key_finished_output(iw, logger, mm_section, &
302 "PRINT%FF_INFO")
303 CALL timestop(handle)
304
305 DEALLOCATE (namearray)
306 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 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 CALL timeset(routinen, handle)
337 NULLIFY (logger)
338 logger => cp_get_default_logger()
339 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
340 extension=".mmLog")
341
342 avail_section(1) = "BOND"; bond_section(1) = avail_section(1)
343 avail_section(11) = "BONDS"
344 avail_section(2) = "ANGL"; angl_section(1) = avail_section(2)
345 avail_section(3) = "THETA"; angl_section(2) = avail_section(3)
346 avail_section(12) = "THETAS"
347 avail_section(13) = "ANGLE"
348 avail_section(14) = "ANGLES"
349 avail_section(4) = "DIHEDRAL"; thet_section(1) = avail_section(4)
350 avail_section(15) = "DIHEDRALS"
351 avail_section(5) = "PHI"; thet_section(2) = avail_section(5)
352 avail_section(6) = "IMPROPER"; impr_section(1) = avail_section(6)
353 avail_section(20) = "IMPROPERS"
354 avail_section(7) = "IMPH"; impr_section(2) = avail_section(7)
355 avail_section(16) = "IMPHI"
356 avail_section(8) = "NONBONDED"; nbon_section(1) = avail_section(8)
357 avail_section(9) = "NBOND"; nbon_section(2) = avail_section(9)
358 avail_section(10) = "HBOND"
359 avail_section(17) = "NBFIX"
360 avail_section(18) = "CMAP"
361 avail_section(19) = "END"
362
363 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 nbond = 0
371 DO ilab = 1, SIZE(bond_section)
372 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
373 label = trim(bond_section(ilab))
374 DO
375 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
376 IF (found) THEN
377 CALL parser_get_object(parser, string)
378 IF (index(string, trim(label)) /= 1) cycle
379 CALL charmm_get_next_line(parser, 1)
380 DO
381 CALL parser_get_object(parser, string)
382 CALL uppercase(string)
383 IF (any(string == avail_section)) EXIT
384 CALL parser_get_object(parser, string2)
385 CALL uppercase(string2)
386 nbond = nbond + 1
387 CALL reallocate(chm_info%bond_a, 1, nbond)
388 CALL reallocate(chm_info%bond_b, 1, nbond)
389 CALL reallocate(chm_info%bond_k, 1, nbond)
390 CALL reallocate(chm_info%bond_r0, 1, nbond)
391 chm_info%bond_a(nbond) = string
392 chm_info%bond_b(nbond) = string2
393 CALL parser_get_object(parser, chm_info%bond_k(nbond))
394 CALL parser_get_object(parser, chm_info%bond_r0(nbond))
395 IF (iw > 0) WRITE (iw, *) " CHM BOND ", nbond, &
396 trim(chm_info%bond_a(nbond)), " ", &
397 trim(chm_info%bond_b(nbond)), " ", &
398 chm_info%bond_k(nbond), &
399 chm_info%bond_r0(nbond)
400 ! Do some units conversion into internal atomic units
401 chm_info%bond_r0(nbond) = cp_unit_to_cp2k(chm_info%bond_r0(nbond), "angstrom")
402 chm_info%bond_k(nbond) = cp_unit_to_cp2k(chm_info%bond_k(nbond), "kcalmol*angstrom^-2")
403 CALL charmm_get_next_line(parser, 1)
404 END DO
405 ELSE
406 EXIT
407 END IF
408 END DO
409 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 nbend = 0
421 nub = 0
422 DO ilab = 1, SIZE(angl_section)
423 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
424 label = trim(angl_section(ilab))
425 DO
426 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
427 IF (found) THEN
428 CALL parser_get_object(parser, string)
429 IF (index(string, trim(label)) /= 1) cycle
430 CALL charmm_get_next_line(parser, 1)
431 DO
432 CALL parser_get_object(parser, string)
433 CALL uppercase(string)
434 IF (any(string == avail_section)) EXIT
435 CALL parser_get_object(parser, string2)
436 CALL parser_get_object(parser, string3)
437 CALL uppercase(string2)
438 CALL uppercase(string3)
439 nbend = nbend + 1
440 CALL reallocate(chm_info%bend_a, 1, nbend)
441 CALL reallocate(chm_info%bend_b, 1, nbend)
442 CALL reallocate(chm_info%bend_c, 1, nbend)
443 CALL reallocate(chm_info%bend_k, 1, nbend)
444 CALL reallocate(chm_info%bend_theta0, 1, nbend)
445 chm_info%bend_a(nbend) = string
446 chm_info%bend_b(nbend) = string2
447 chm_info%bend_c(nbend) = string3
448 CALL parser_get_object(parser, chm_info%bend_k(nbend))
449 CALL parser_get_object(parser, chm_info%bend_theta0(nbend))
450 IF (iw > 0) WRITE (iw, *) " CHM BEND ", nbend, &
451 trim(chm_info%bend_a(nbend)), " ", &
452 trim(chm_info%bend_b(nbend)), " ", &
453 trim(chm_info%bend_c(nbend)), " ", &
454 chm_info%bend_k(nbend), &
455 chm_info%bend_theta0(nbend)
456 ! Do some units conversion into internal atomic units
457 chm_info%bend_theta0(nbend) = cp_unit_to_cp2k(chm_info%bend_theta0(nbend), "deg")
458 chm_info%bend_k(nbend) = cp_unit_to_cp2k(chm_info%bend_k(nbend), "kcalmol*rad^-2")
459 IF (parser_test_next_token(parser) == "FLT") THEN
460 nub = nub + 1
461 CALL reallocate(chm_info%ub_a, 1, nub)
462 CALL reallocate(chm_info%ub_b, 1, nub)
463 CALL reallocate(chm_info%ub_c, 1, nub)
464 CALL reallocate(chm_info%ub_k, 1, nub)
465 CALL reallocate(chm_info%ub_r0, 1, nub)
466 chm_info%ub_a(nub) = string
467 chm_info%ub_b(nub) = string2
468 chm_info%ub_c(nub) = string3
469 CALL parser_get_object(parser, chm_info%ub_k(nub))
470 CALL parser_get_object(parser, chm_info%ub_r0(nub))
471 IF (iw > 0) WRITE (iw, *) " CHM UB ", nub, &
472 trim(chm_info%ub_a(nub)), " ", &
473 trim(chm_info%ub_b(nub)), " ", &
474 trim(chm_info%ub_c(nub)), " ", &
475 chm_info%ub_k(nub), &
476 chm_info%ub_r0(nub)
477 ! Do some units conversion into internal atomic units
478 chm_info%ub_r0(nub) = cp_unit_to_cp2k(chm_info%ub_r0(nub), "angstrom")
479 chm_info%ub_k(nub) = cp_unit_to_cp2k(chm_info%ub_k(nub), "kcalmol*angstrom^-2")
480 END IF
481 CALL charmm_get_next_line(parser, 1)
482 END DO
483 ELSE
484 EXIT
485 END IF
486 END DO
487 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 ntorsion = 0
496 DO ilab = 1, SIZE(thet_section)
497 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
498 label = trim(thet_section(ilab))
499 DO
500 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
501 IF (found) THEN
502 CALL parser_get_object(parser, string)
503 IF (index(string, trim(label)) /= 1) cycle
504 CALL charmm_get_next_line(parser, 1)
505 DO
506 CALL parser_get_object(parser, string)
507 CALL uppercase(string)
508 IF (any(string == avail_section)) EXIT
509 CALL parser_get_object(parser, string2)
510 CALL parser_get_object(parser, string3)
511 CALL parser_get_object(parser, string4)
512 CALL uppercase(string2)
513 CALL uppercase(string3)
514 CALL uppercase(string4)
515 ntorsion = ntorsion + 1
516 CALL reallocate(chm_info%torsion_a, 1, ntorsion)
517 CALL reallocate(chm_info%torsion_b, 1, ntorsion)
518 CALL reallocate(chm_info%torsion_c, 1, ntorsion)
519 CALL reallocate(chm_info%torsion_d, 1, ntorsion)
520 CALL reallocate(chm_info%torsion_k, 1, ntorsion)
521 CALL reallocate(chm_info%torsion_m, 1, ntorsion)
522 CALL reallocate(chm_info%torsion_phi0, 1, ntorsion)
523 chm_info%torsion_a(ntorsion) = string
524 chm_info%torsion_b(ntorsion) = string2
525 chm_info%torsion_c(ntorsion) = string3
526 chm_info%torsion_d(ntorsion) = string4
527 CALL parser_get_object(parser, chm_info%torsion_k(ntorsion))
528 CALL parser_get_object(parser, chm_info%torsion_m(ntorsion))
529 CALL parser_get_object(parser, chm_info%torsion_phi0(ntorsion))
530 IF (iw > 0) WRITE (iw, *) " CHM TORSION ", ntorsion, &
531 trim(chm_info%torsion_a(ntorsion)), " ", &
532 trim(chm_info%torsion_b(ntorsion)), " ", &
533 trim(chm_info%torsion_c(ntorsion)), " ", &
534 trim(chm_info%torsion_d(ntorsion)), " ", &
535 chm_info%torsion_k(ntorsion), &
536 chm_info%torsion_m(ntorsion), &
537 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 "deg")
541 chm_info%torsion_k(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_k(ntorsion), "kcalmol")
542 CALL charmm_get_next_line(parser, 1)
543 END DO
544 ELSE
545 EXIT
546 END IF
547 END DO
548 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 nimpr = 0
557 DO ilab = 1, SIZE(impr_section)
558 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
559 label = trim(impr_section(ilab))
560 DO
561 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
562 IF (found) THEN
563 CALL parser_get_object(parser, string)
564 IF (index(string, trim(label)) /= 1) cycle
565 CALL charmm_get_next_line(parser, 1)
566 DO
567 CALL parser_get_object(parser, string)
568 CALL uppercase(string)
569 IF (any(string == avail_section)) EXIT
570 CALL parser_get_object(parser, string2)
571 CALL parser_get_object(parser, string3)
572 CALL parser_get_object(parser, string4)
573 CALL uppercase(string2)
574 CALL uppercase(string3)
575 CALL uppercase(string4)
576 nimpr = nimpr + 1
577 CALL reallocate(chm_info%impr_a, 1, nimpr)
578 CALL reallocate(chm_info%impr_b, 1, nimpr)
579 CALL reallocate(chm_info%impr_c, 1, nimpr)
580 CALL reallocate(chm_info%impr_d, 1, nimpr)
581 CALL reallocate(chm_info%impr_k, 1, nimpr)
582 CALL reallocate(chm_info%impr_phi0, 1, nimpr)
583 chm_info%impr_a(nimpr) = string
584 chm_info%impr_b(nimpr) = string2
585 chm_info%impr_c(nimpr) = string3
586 chm_info%impr_d(nimpr) = string4
587 CALL parser_get_object(parser, chm_info%impr_k(nimpr))
588 CALL parser_get_object(parser, dummy)
589 CALL parser_get_object(parser, chm_info%impr_phi0(nimpr))
590 IF (iw > 0) WRITE (iw, *) " CHM IMPROPERS ", nimpr, &
591 trim(chm_info%impr_a(nimpr)), " ", &
592 trim(chm_info%impr_b(nimpr)), " ", &
593 trim(chm_info%impr_c(nimpr)), " ", &
594 trim(chm_info%impr_d(nimpr)), " ", &
595 chm_info%impr_k(nimpr), &
596 chm_info%impr_phi0(nimpr)
597 ! Do some units conversion into internal atomic units
598 chm_info%impr_phi0(nimpr) = cp_unit_to_cp2k(chm_info%impr_phi0(nimpr), "deg")
599 chm_info%impr_k(nimpr) = cp_unit_to_cp2k(chm_info%impr_k(nimpr), "kcalmol")
600 CALL charmm_get_next_line(parser, 1)
601 END DO
602 ELSE
603 EXIT
604 END IF
605 END DO
606 CALL parser_release(parser)
607 END DO
608 !-----------------------------------------------------------------------------
609 !-----------------------------------------------------------------------------
610 ! 5. Read in all the Nonbonded info from the param file here
611 !-----------------------------------------------------------------------------
612 nnonbond = 0
613 nonfo = 0
614 DO ilab = 1, SIZE(nbon_section)
615 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env)
616 label = trim(nbon_section(ilab))
617 DO
618 CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
619 IF (found) THEN
620 CALL parser_get_object(parser, string)
621 IF (index(string, trim(label)) /= 1) cycle
622 CALL charmm_get_next_line(parser, 1)
623 DO
624 CALL parser_get_object(parser, string)
625 CALL uppercase(string)
626 IF (any(string == avail_section)) EXIT
627 nnonbond = nnonbond + 1
628 CALL reallocate(chm_info%nonbond_a, 1, nnonbond)
629 CALL reallocate(chm_info%nonbond_eps, 1, nnonbond)
630 CALL reallocate(chm_info%nonbond_rmin2, 1, nnonbond)
631 chm_info%nonbond_a(nnonbond) = string
632 CALL parser_get_object(parser, chm_info%nonbond_eps(nnonbond))
633 CALL parser_get_object(parser, chm_info%nonbond_eps(nnonbond))
634 CALL parser_get_object(parser, chm_info%nonbond_rmin2(nnonbond))
635 IF (iw > 0) WRITE (iw, *) " CHM NONBOND ", nnonbond, &
636 trim(chm_info%nonbond_a(nnonbond)), " ", &
637 chm_info%nonbond_eps(nnonbond), &
638 chm_info%nonbond_rmin2(nnonbond)
639 chm_info%nonbond_rmin2(nnonbond) = cp_unit_to_cp2k(chm_info%nonbond_rmin2(nnonbond), &
640 "angstrom")
641 chm_info%nonbond_eps(nnonbond) = cp_unit_to_cp2k(chm_info%nonbond_eps(nnonbond), &
642 "kcalmol")
643 IF (parser_test_next_token(parser) == "FLT") THEN
644 nonfo = nonfo + 1
645 CALL reallocate(chm_info%nonbond_a_14, 1, nonfo)
646 CALL reallocate(chm_info%nonbond_eps_14, 1, nonfo)
647 CALL reallocate(chm_info%nonbond_rmin2_14, 1, nonfo)
648 chm_info%nonbond_a_14(nonfo) = chm_info%nonbond_a(nnonbond)
649 CALL parser_get_object(parser, chm_info%nonbond_eps_14(nonfo))
650 CALL parser_get_object(parser, chm_info%nonbond_eps_14(nonfo))
651 CALL parser_get_object(parser, chm_info%nonbond_rmin2_14(nonfo))
652 IF (iw > 0) WRITE (iw, *) " CHM ONFO ", nonfo, &
653 trim(chm_info%nonbond_a_14(nonfo)), " ", &
654 chm_info%nonbond_eps_14(nonfo), &
655 chm_info%nonbond_rmin2_14(nonfo)
656 chm_info%nonbond_rmin2_14(nonfo) = cp_unit_to_cp2k(chm_info%nonbond_rmin2_14(nonfo), &
657 "angstrom")
658 chm_info%nonbond_eps_14(nonfo) = cp_unit_to_cp2k(chm_info%nonbond_eps_14(nonfo), &
659 "kcalmol")
660 END IF
661 CALL charmm_get_next_line(parser, 1)
662 END DO
663 ELSE
664 EXIT
665 END IF
666 END DO
667 CALL parser_release(parser)
668 END DO
669 CALL cp_print_key_finished_output(iw, logger, mm_section, &
670 "PRINT%FF_INFO")
671 CALL timestop(handle)
672
673 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 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 CALL timeset(routinen, handle)
697 NULLIFY (logger)
698 logger => cp_get_default_logger()
699 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
700 extension=".mmLog")
701
702 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 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 DO i = 1, SIZE(amb_info%bond_a)
714 IF (iw > 0) WRITE (iw, *) " AMB BOND ", i, &
715 trim(amb_info%bond_a(i)), " ", &
716 trim(amb_info%bond_b(i)), " ", &
717 amb_info%bond_k(i), &
718 amb_info%bond_r0(i)
719
720 ! Do some units conversion into internal atomic units
721 amb_info%bond_r0(i) = cp_unit_to_cp2k(amb_info%bond_r0(i), "angstrom")
722 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 DO i = 1, SIZE(amb_info%bend_a)
734 IF (iw > 0) WRITE (iw, *) " AMB BEND ", i, &
735 trim(amb_info%bend_a(i)), " ", &
736 trim(amb_info%bend_b(i)), " ", &
737 trim(amb_info%bend_c(i)), " ", &
738 amb_info%bend_k(i), &
739 amb_info%bend_theta0(i)
740
741 ! Do some units conversion into internal atomic units
742 amb_info%bend_theta0(i) = cp_unit_to_cp2k(amb_info%bend_theta0(i), "rad")
743 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 DO i = 1, SIZE(amb_info%torsion_a)
752 IF (iw > 0) WRITE (iw, *) " AMB TORSION ", i, &
753 trim(amb_info%torsion_a(i)), " ", &
754 trim(amb_info%torsion_b(i)), " ", &
755 trim(amb_info%torsion_c(i)), " ", &
756 trim(amb_info%torsion_d(i)), " ", &
757 amb_info%torsion_k(i), &
758 amb_info%torsion_m(i), &
759 amb_info%torsion_phi0(i)
760
761 ! Do some units conversion into internal atomic units
762 amb_info%torsion_phi0(i) = cp_unit_to_cp2k(amb_info%torsion_phi0(i), "rad")
763 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 IF (ASSOCIATED(amb_info%raw_torsion_phi0)) THEN
768 DO i = 1, SIZE(amb_info%raw_torsion_k)
769 amb_info%raw_torsion_phi0(i) = cp_unit_to_cp2k(amb_info%raw_torsion_phi0(i), "rad")
770 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 DO i = 1, SIZE(amb_info%nonbond_eps)
778 IF (iw > 0) WRITE (iw, *) " AMB NONBOND ", i, &
779 trim(amb_info%nonbond_a(i)), " ", &
780 amb_info%nonbond_eps(i), &
781 amb_info%nonbond_rmin2(i)
782
783 ! Do some units conversion into internal atomic units
784 amb_info%nonbond_rmin2(i) = cp_unit_to_cp2k(amb_info%nonbond_rmin2(i), "angstrom")
785 amb_info%nonbond_eps(i) = cp_unit_to_cp2k(amb_info%nonbond_eps(i), "kcalmol")
786 END DO
787 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
788 CALL timestop(handle)
789 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 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 DO i = 1, nline
810 len_line = len_trim(parser%input_line)
811 DO WHILE (parser%input_line(len_line:len_line) == continuation_char)
812 CALL parser_get_next_line(parser, 1)
813 len_line = len_trim(parser%input_line)
814 END DO
815 CALL parser_get_next_line(parser, 1)
816 END DO
817
818 END SUBROUTINE charmm_get_next_line
819
820END MODULE force_fields_ext
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
character(len=3) function, public parser_test_next_token(parser, string_length)
Test next input object.
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1150
Define all structure types related to force field kinds.
integer, parameter, public do_ff_g96
Define all structures types related to force_fields.
subroutine, public read_force_field_gromos(ff_type, para_env, mm_section)
Reads the GROMOS force_field.
subroutine, public read_force_field_charmm(ff_type, para_env, mm_section)
Reads the charmm force_field.
subroutine, public read_force_field_amber(ff_type, para_env, mm_section, particle_set)
Reads the AMBER force_field.
objects that represent the structure of input sections and the data contained in an input section
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Utility routines for the memory handling.
Interface to the message passing library MPI.
Define the data structure for the particle information.
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Handles all functions used to read and interpret AMBER coordinates and topology files.
subroutine, public rdparm_amber_8(filename, output_unit, para_env, do_connectivity, do_forcefield, atom_info, conn_info, amb_info, particle_set)
Access information form the AMBER topology file Notes on file structure:
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment