(git:34ef472)
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 ! **************************************************************************************************
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  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 
820 END MODULE force_fields_ext
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
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)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
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,...
Definition: cp_output_handling.F:1013
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
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 ...
Definition: cp_parser_methods.F:421
cp_parser_methods::parser_test_next_token
character(len=3) function, public parser_test_next_token(parser, string_length)
Test next input object.
Definition: cp_parser_methods.F:691
cp_parser_methods::parser_search_string
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 ...
Definition: cp_parser_methods.F:786
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition: cp_parser_types.F:109
cp_parser_types::parser_create
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.
Definition: cp_parser_types.F:147
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_to_cp2k
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
force_field_kind_types
Define all structure types related to force field kinds.
Definition: force_field_kind_types.F:14
force_field_kind_types::do_ff_g96
integer, parameter, public do_ff_g96
Definition: force_field_kind_types.F:25
force_field_types
Define all structures types related to force_fields.
Definition: force_field_types.F:16
force_fields_ext
Definition: force_fields_ext.F:22
force_fields_ext::read_force_field_gromos
subroutine, public read_force_field_gromos(ff_type, para_env, mm_section)
Reads the GROMOS force_field.
Definition: force_fields_ext.F:69
force_fields_ext::read_force_field_charmm
subroutine, public read_force_field_charmm(ff_type, para_env, mm_section)
Reads the charmm force_field.
Definition: force_fields_ext.F:316
force_fields_ext::read_force_field_amber
subroutine, public read_force_field_amber(ff_type, para_env, mm_section, particle_set)
Reads the AMBER force_field.
Definition: force_fields_ext.F:684
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition: string_utilities.F:3361
topology_amber
Handles all functions used to read and interpret AMBER coordinates and topology files.
Definition: topology_amber.F:14
topology_amber::rdparm_amber_8
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:
Definition: topology_amber.F:368