(git:1f285aa)
force_fields.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 !> \author CJM
18 ! **************************************************************************************************
19 MODULE force_fields
20  USE atomic_kind_types, ONLY: atomic_kind_type
21  USE cell_types, ONLY: cell_type
22  USE cp_log_handling, ONLY: cp_get_default_logger,&
23  cp_logger_type
24  USE cp_output_handling, ONLY: cp_p_file,&
25  cp_print_key_finished_output,&
26  cp_print_key_should_output,&
27  cp_print_key_unit_nr
28  USE cp_units, ONLY: cp_unit_from_cp2k
29  USE ewald_environment_types, ONLY: ewald_environment_type
30  USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type
31  USE force_field_kind_types, ONLY: do_ff_amber,&
32  do_ff_charmm,&
33  do_ff_g87,&
34  do_ff_g96,&
35  do_ff_undef
36  USE force_field_types, ONLY: deallocate_ff_type,&
37  force_field_type,&
38  init_ff_type
39  USE force_fields_ext, ONLY: read_force_field_amber,&
40  read_force_field_charmm,&
41  read_force_field_gromos
42  USE force_fields_input, ONLY: read_force_field_section
43  USE force_fields_util, ONLY: clean_intra_force_kind,&
44  force_field_pack,&
45  force_field_qeff_output
46  USE input_constants, ONLY: do_skip_14
47  USE input_section_types, ONLY: section_vals_get_subs_vals,&
48  section_vals_type,&
49  section_vals_val_get
50  USE kinds, ONLY: dp
51  USE message_passing, ONLY: mp_para_env_type
52  USE molecule_kind_types, ONLY: molecule_kind_type
53  USE molecule_types, ONLY: molecule_type
54  USE particle_types, ONLY: particle_type
55  USE qmmm_types_low, ONLY: qmmm_env_mm_type
56 #include "./base/base_uses.f90"
57 
58  IMPLICIT NONE
59 
60  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields'
61 
62  PRIVATE
63  PUBLIC :: force_field_control
64 
65 CONTAINS
66 
67 ! **************************************************************************************************
68 !> \brief 1. If reading in from external file, make sure its there first
69 !> 2. Read in the force_field from the corresponding locations
70 !> \param atomic_kind_set ...
71 !> \param particle_set ...
72 !> \param molecule_kind_set ...
73 !> \param molecule_set ...
74 !> \param ewald_env ...
75 !> \param fist_nonbond_env ...
76 !> \param root_section ...
77 !> \param para_env ...
78 !> \param qmmm ...
79 !> \param qmmm_env ...
80 !> \param subsys_section ...
81 !> \param mm_section ...
82 !> \param shell_particle_set ...
83 !> \param core_particle_set ...
84 !> \param cell ...
85 ! **************************************************************************************************
86  SUBROUTINE force_field_control(atomic_kind_set, particle_set, &
87  molecule_kind_set, molecule_set, ewald_env, fist_nonbond_env, &
88  root_section, para_env, qmmm, qmmm_env, subsys_section, mm_section, &
89  shell_particle_set, core_particle_set, cell)
90 
91  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
92  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
93  TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
94  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
95  TYPE(ewald_environment_type), POINTER :: ewald_env
96  TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
97  TYPE(section_vals_type), POINTER :: root_section
98  TYPE(mp_para_env_type), POINTER :: para_env
99  LOGICAL, INTENT(IN), OPTIONAL :: qmmm
100  TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
101  TYPE(section_vals_type), POINTER :: subsys_section, mm_section
102  TYPE(particle_type), DIMENSION(:), POINTER :: shell_particle_set, core_particle_set
103  TYPE(cell_type), POINTER :: cell
104 
105  CHARACTER(len=*), PARAMETER :: routinen = 'force_field_control'
106 
107  INTEGER :: exclude_ei, exclude_vdw, handle, iw
108  LOGICAL :: found
109  TYPE(cp_logger_type), POINTER :: logger
110  TYPE(force_field_type) :: ff_type
111  TYPE(section_vals_type), POINTER :: topology_section
112 
113  CALL timeset(routinen, handle)
114  logger => cp_get_default_logger()
115 
116  iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
117  extension=".mmLog")
118 
119  !-----------------------------------------------------------------------------
120  ! 1. Initialize the ff_type structure type
121  !-----------------------------------------------------------------------------
122  CALL init_ff_type(ff_type)
123 
124  !-----------------------------------------------------------------------------
125  ! 2. Read in the force field section in the input file if any
126  !-----------------------------------------------------------------------------
127  CALL read_force_field_section(ff_type, para_env, mm_section)
128 
129  !-----------------------------------------------------------------------------
130  ! 2.1 In case exclusion 1-4 was requested, we need to modify the values of
131  ! the scale factors setting them to zero..
132  !-----------------------------------------------------------------------------
133  topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
134  CALL section_vals_val_get(topology_section, "EXCLUDE_VDW", i_val=exclude_vdw)
135  CALL section_vals_val_get(topology_section, "EXCLUDE_EI", i_val=exclude_ei)
136  IF (exclude_vdw == do_skip_14) ff_type%vdw_scale14 = 0.0_dp
137  IF (exclude_ei == do_skip_14) ff_type%ei_scale14 = 0.0_dp
138 
139  !-----------------------------------------------------------------------------
140  ! 3. If reading in from external file, make sure its there first
141  !-----------------------------------------------------------------------------
142  SELECT CASE (ff_type%ff_type)
143  CASE (do_ff_charmm, do_ff_amber, do_ff_g96, do_ff_g87)
144  INQUIRE (file=ff_type%ff_file_name, exist=found)
145  IF (.NOT. found) THEN
146  cpabort("Force field file missing")
147  END IF
148  CASE (do_ff_undef)
149  ! Do Nothing
150  CASE DEFAULT
151  cpabort("Force field type not implemented")
152  END SELECT
153 
154  !-----------------------------------------------------------------------------
155  ! 4. Read in the force field from the corresponding locations
156  !-----------------------------------------------------------------------------
157  SELECT CASE (ff_type%ff_type)
158  CASE (do_ff_charmm)
159  CALL read_force_field_charmm(ff_type, para_env, mm_section)
160  CASE (do_ff_amber)
161  CALL read_force_field_amber(ff_type, para_env, mm_section, particle_set)
162  CASE (do_ff_g87, do_ff_g96)
163  CALL read_force_field_gromos(ff_type, para_env, mm_section)
164  CASE (do_ff_undef)
165  ! Do Nothing
166  CASE DEFAULT
167  cpabort("Force field type not implemented")
168  END SELECT
169 
170  !-----------------------------------------------------------------------------
171  ! 5. Possibly print the top file
172  !-----------------------------------------------------------------------------
173  CALL print_pot_parameter_file(ff_type, mm_section)
174 
175  !-----------------------------------------------------------------------------
176  ! 6. Pack all force field info into different structures
177  !-----------------------------------------------------------------------------
178  CALL force_field_pack(particle_set, atomic_kind_set, molecule_kind_set, molecule_set, &
179  ewald_env, fist_nonbond_env, ff_type, root_section, qmmm, qmmm_env, mm_section, &
180  subsys_section, shell_particle_set=shell_particle_set, &
181  core_particle_set=core_particle_set, cell=cell)
182 
183  !-----------------------------------------------------------------------------
184  ! 7. Output total system charge assigned to qeff
185  !-----------------------------------------------------------------------------
186  CALL force_field_qeff_output(particle_set, molecule_kind_set, &
187  molecule_set, mm_section, fist_nonbond_env%charges)
188 
189  !-----------------------------------------------------------------------------
190  ! 8. Clean up "UNSET" bond,bend,UB,TORSION,IMPR,ONFO kinds
191  !-----------------------------------------------------------------------------
192  CALL clean_intra_force_kind(molecule_kind_set, mm_section)
193 
194  !-----------------------------------------------------------------------------
195  ! 9. Cleanup the ff_type structure type
196  !-----------------------------------------------------------------------------
197  CALL deallocate_ff_type(ff_type)
198 
199  CALL cp_print_key_finished_output(iw, logger, mm_section, &
200  "PRINT%FF_INFO")
201  CALL timestop(handle)
202 
203  END SUBROUTINE force_field_control
204 
205 ! **************************************************************************************************
206 !> \brief Prints force field information in a pot file
207 !> \param ff_type ...
208 !> \param mm_section ...
209 !> \author Teodoro Laino [tlaino, teodoro.laino-AT-gmail.com] - 11.2008
210 ! **************************************************************************************************
211  SUBROUTINE print_pot_parameter_file(ff_type, mm_section)
212 
213  TYPE(force_field_type) :: ff_type
214  TYPE(section_vals_type), POINTER :: mm_section
215 
216  CHARACTER(len=*), PARAMETER :: routinen = 'print_pot_parameter_file'
217 
218  INTEGER :: handle, i, iw, m
219  REAL(kind=dp) :: eps, k, phi0, r0, sigma, theta0
220  TYPE(cp_logger_type), POINTER :: logger
221 
222  CALL timeset(routinen, handle)
223  logger => cp_get_default_logger()
224  IF (btest(cp_print_key_should_output(logger%iter_info, mm_section, "PRINT%FF_PARAMETER_FILE") &
225  , cp_p_file)) THEN
226  iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_PARAMETER_FILE", &
227  middle_name="force_field", extension=".pot")
228  IF (iw > 0) THEN
229  ! Header
230  WRITE (iw, 1000) "Force Field Parameter File dumped into CHARMM FF style"
231  END IF
232  SELECT CASE (ff_type%ff_type)
233  CASE (do_ff_charmm)
234  cpwarn("Dumping FF parameter file for CHARMM FF not implemented!")
235  CASE (do_ff_amber)
236  IF (iw > 0) THEN
237  ! Bonds
238  WRITE (iw, 1001)
239  DO i = 1, SIZE(ff_type%amb_info%bond_a)
240  k = cp_unit_from_cp2k(ff_type%amb_info%bond_k(i), "kcalmol*angstrom^-2")
241  r0 = cp_unit_from_cp2k(ff_type%amb_info%bond_r0(i), "angstrom")
242  WRITE (iw, 2001) ff_type%amb_info%bond_a(i), &
243  ff_type%amb_info%bond_b(i), &
244  k, r0
245  END DO
246  ! Angles
247  WRITE (iw, 1002)
248  DO i = 1, SIZE(ff_type%amb_info%bend_a)
249  k = cp_unit_from_cp2k(ff_type%amb_info%bend_k(i), "kcalmol*rad^-2")
250  theta0 = cp_unit_from_cp2k(ff_type%amb_info%bend_theta0(i), "deg")
251  WRITE (iw, 2002) ff_type%amb_info%bend_a(i), &
252  ff_type%amb_info%bend_b(i), &
253  ff_type%amb_info%bend_c(i), &
254  k, theta0
255  END DO
256  ! Torsions
257  WRITE (iw, 1003)
258  DO i = 1, SIZE(ff_type%amb_info%torsion_a)
259  k = cp_unit_from_cp2k(ff_type%amb_info%torsion_k(i), "kcalmol")
260  m = ff_type%amb_info%torsion_m(i)
261  phi0 = cp_unit_from_cp2k(ff_type%amb_info%torsion_phi0(i), "deg")
262  WRITE (iw, 2003) ff_type%amb_info%torsion_a(i), &
263  ff_type%amb_info%torsion_b(i), &
264  ff_type%amb_info%torsion_c(i), &
265  ff_type%amb_info%torsion_d(i), &
266  k, m, phi0
267  END DO
268  ! Lennard-Jones
269  WRITE (iw, 1005)
270  DO i = 1, SIZE(ff_type%amb_info%nonbond_a)
271  eps = cp_unit_from_cp2k(ff_type%amb_info%nonbond_eps(i), "kcalmol")
272  sigma = cp_unit_from_cp2k(ff_type%amb_info%nonbond_rmin2(i), "angstrom")
273  WRITE (iw, 2005) ff_type%amb_info%nonbond_a(i), &
274  eps, sigma
275  END DO
276  END IF
277  CASE (do_ff_g87, do_ff_g96)
278  cpwarn("Dumping FF parameter file for GROMOS FF not implemented!")
279  CASE (do_ff_undef)
280  cpwarn("Dumping FF parameter file for INPUT FF not implemented!")
281  END SELECT
282  IF (iw > 0) THEN
283  WRITE (iw, '(/,A)') "END"
284  END IF
285  CALL cp_print_key_finished_output(iw, logger, mm_section, &
286  "PRINT%FF_PARAMETER_FILE")
287  END IF
288  CALL timestop(handle)
289  RETURN
290 1000 FORMAT("*>>>>>>>", t12, a, t73, "<<<<<<<")
291 1001 FORMAT(/, "BONDS", /, "!", /, "!V(bond) = Kb(b - b0)**2", /, "!", /, "!Kb: kcal/mole/A**2", /, &
292  "!b0: A", /, "!", /, "! atom type Kb b0", /, "!")
293 1002 FORMAT(/, "ANGLES", /, "!", /, "!V(angle) = Ktheta(Theta - Theta0)**2", /, "!", /, &
294  "!V(Urey-Bradley) = Kub(S - S0)**2", /, "!", /, "!Ktheta: kcal/mole/rad**2", /, &
295  "!Theta0: degrees", /, "!Kub: kcal/mole/A**2 (Urey-Bradley)", /, "!S0: A", /, &
296  "!", /, "! atom types Ktheta Theta0 Kub S0", /, "!")
297 1003 FORMAT(/, "DIHEDRALS", /, "!", /, "!V(dihedral) = Kchi(1 + cos(n(chi) - delta))", /, &
298  "!", /, "!Kchi: kcal/mole", /, "!n: multiplicity", /, "!delta: degrees", /, &
299  "!", /, "! atom types Kchi n delta", /, "!")
300 1005 FORMAT(/, "NONBONDED", /, "!", /, &
301  "!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]", /, &
302  "!", /, "!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)", /, &
303  "!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j", /, "!", /, &
304  "!atom ignored epsilon Rmin/2 ignored eps,1-4 "// &
305  "Rmin/2,1-4", /, "!")
306 
307 2001 FORMAT(a6, 1x, a6, 1x, 2f15.9) ! bond
308 2002 FORMAT(a6, 1x, a6, 1x, a6, 1x, 2f15.9) ! angle
309 2003 FORMAT(a6, 1x, a6, 1x, a6, 1x, a6, 1x, f15.9, i5, f15.9) ! torsion
310 2005 FORMAT(a6, 1x, " 0.000000000", 2f15.9) ! nonbond
311  END SUBROUTINE print_pot_parameter_file
312 
313 END MODULE force_fields
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
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_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
ewald_environment_types
Definition: ewald_environment_types.F:14
fist_nonbond_env_types
Definition: fist_nonbond_env_types.F:13
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_undef
integer, parameter, public do_ff_undef
Definition: force_field_kind_types.F:25
force_field_kind_types::do_ff_charmm
integer, parameter, public do_ff_charmm
Definition: force_field_kind_types.F:25
force_field_kind_types::do_ff_g87
integer, parameter, public do_ff_g87
Definition: force_field_kind_types.F:25
force_field_kind_types::do_ff_g96
integer, parameter, public do_ff_g96
Definition: force_field_kind_types.F:25
force_field_kind_types::do_ff_amber
integer, parameter, public do_ff_amber
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_field_types::init_ff_type
subroutine, public init_ff_type(ff_type)
Just NULLIFY and zero all the stuff
Definition: force_field_types.F:231
force_field_types::deallocate_ff_type
subroutine, public deallocate_ff_type(ff_type)
Just DEALLOCATE all the stuff
Definition: force_field_types.F:450
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
force_fields_input
Definition: force_fields_input.F:20
force_fields_input::read_force_field_section
subroutine, public read_force_field_section(ff_type, para_env, mm_section)
Reads the force_field input section.
Definition: force_fields_input.F:2345
force_fields_util
Definition: force_fields_util.F:19
force_fields_util::force_field_pack
subroutine, public force_field_pack(particle_set, atomic_kind_set, molecule_kind_set, molecule_set, ewald_env, fist_nonbond_env, ff_type, root_section, qmmm, qmmm_env, mm_section, subsys_section, shell_particle_set, core_particle_set, cell)
Pack in all the information needed for the force_field.
Definition: force_fields_util.F:111
force_fields_util::force_field_qeff_output
subroutine, public force_field_qeff_output(particle_set, molecule_kind_set, molecule_set, mm_section, charges)
Compute the total qeff charges for each molecule kind and total system.
Definition: force_fields_util.F:401
force_fields_util::clean_intra_force_kind
subroutine, public clean_intra_force_kind(molecule_kind_set, mm_section)
Removes UNSET force field types.
Definition: force_fields_util.F:496
force_fields
Definition: force_fields.F:19
force_fields::force_field_control
subroutine, public force_field_control(atomic_kind_set, particle_set, molecule_kind_set, molecule_set, ewald_env, fist_nonbond_env, root_section, para_env, qmmm, qmmm_env, subsys_section, mm_section, shell_particle_set, core_particle_set, cell)
If reading in from external file, make sure its there first
Definition: force_fields.F:90
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_skip_14
integer, parameter, public do_skip_14
Definition: input_constants.F:412
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
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition: molecule_kind_types.F:16
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qmmm_types_low
Definition: qmmm_types_low.F:13