Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 7929 : 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 2643 : CALL timeset(routineN, handle)
114 2643 : logger => cp_get_default_logger()
115 :
116 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
117 2643 : extension=".mmLog")
118 :
119 : !-----------------------------------------------------------------------------
120 : ! 1. Initialize the ff_type structure type
121 : !-----------------------------------------------------------------------------
122 2643 : CALL init_ff_type(ff_type)
123 :
124 : !-----------------------------------------------------------------------------
125 : ! 2. Read in the force field section in the input file if any
126 : !-----------------------------------------------------------------------------
127 2643 : 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 2643 : topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
134 2643 : CALL section_vals_val_get(topology_section, "EXCLUDE_VDW", i_val=exclude_vdw)
135 2643 : CALL section_vals_val_get(topology_section, "EXCLUDE_EI", i_val=exclude_ei)
136 2643 : IF (exclude_vdw == do_skip_14) ff_type%vdw_scale14 = 0.0_dp
137 2643 : 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 3549 : SELECT CASE (ff_type%ff_type)
143 : CASE (do_ff_charmm, do_ff_amber, do_ff_g96, do_ff_g87)
144 906 : INQUIRE (FILE=ff_type%ff_file_name, EXIST=found)
145 906 : IF (.NOT. found) THEN
146 0 : CPABORT("Force field file missing")
147 : END IF
148 : CASE (do_ff_undef)
149 : ! Do Nothing
150 : CASE DEFAULT
151 2643 : CPABORT("Force field type not implemented")
152 : END SELECT
153 :
154 : !-----------------------------------------------------------------------------
155 : ! 4. Read in the force field from the corresponding locations
156 : !-----------------------------------------------------------------------------
157 3523 : SELECT CASE (ff_type%ff_type)
158 : CASE (do_ff_charmm)
159 880 : CALL read_force_field_charmm(ff_type, para_env, mm_section)
160 : CASE (do_ff_amber)
161 14 : CALL read_force_field_amber(ff_type, para_env, mm_section, particle_set)
162 : CASE (do_ff_g87, do_ff_g96)
163 12 : CALL read_force_field_gromos(ff_type, para_env, mm_section)
164 : CASE (do_ff_undef)
165 : ! Do Nothing
166 : CASE DEFAULT
167 2643 : CPABORT("Force field type not implemented")
168 : END SELECT
169 :
170 : !-----------------------------------------------------------------------------
171 : ! 5. Possibly print the top file
172 : !-----------------------------------------------------------------------------
173 2643 : 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 2643 : 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 2643 : molecule_set, mm_section, fist_nonbond_env%charges)
188 :
189 : !-----------------------------------------------------------------------------
190 : ! 8. Clean up "UNSET" bond,bend,UB,TORSION,IMPR,ONFO kinds
191 : !-----------------------------------------------------------------------------
192 2643 : CALL clean_intra_force_kind(molecule_kind_set, mm_section)
193 :
194 : !-----------------------------------------------------------------------------
195 : ! 9. Cleanup the ff_type structure type
196 : !-----------------------------------------------------------------------------
197 2643 : CALL deallocate_ff_type(ff_type)
198 :
199 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
200 2643 : "PRINT%FF_INFO")
201 2643 : CALL timestop(handle)
202 :
203 2643 : 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 2643 : 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 2643 : CALL timeset(routineN, handle)
223 2643 : logger => cp_get_default_logger()
224 2643 : 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 2 : middle_name="force_field", extension=".pot")
228 2 : IF (iw > 0) THEN
229 : ! Header
230 1 : WRITE (iw, 1000) "Force Field Parameter File dumped into CHARMM FF style"
231 : END IF
232 2 : SELECT CASE (ff_type%ff_type)
233 : CASE (do_ff_charmm)
234 0 : CPWARN("Dumping FF parameter file for CHARMM FF not implemented!")
235 : CASE (do_ff_amber)
236 2 : IF (iw > 0) THEN
237 : ! Bonds
238 1 : WRITE (iw, 1001)
239 14 : DO i = 1, SIZE(ff_type%amb_info%bond_a)
240 13 : k = cp_unit_from_cp2k(ff_type%amb_info%bond_k(i), "kcalmol*angstrom^-2")
241 13 : r0 = cp_unit_from_cp2k(ff_type%amb_info%bond_r0(i), "angstrom")
242 13 : WRITE (iw, 2001) ff_type%amb_info%bond_a(i), &
243 13 : ff_type%amb_info%bond_b(i), &
244 27 : k, r0
245 : END DO
246 : ! Angles
247 1 : WRITE (iw, 1002)
248 23 : DO i = 1, SIZE(ff_type%amb_info%bend_a)
249 22 : k = cp_unit_from_cp2k(ff_type%amb_info%bend_k(i), "kcalmol*rad^-2")
250 22 : theta0 = cp_unit_from_cp2k(ff_type%amb_info%bend_theta0(i), "deg")
251 22 : WRITE (iw, 2002) ff_type%amb_info%bend_a(i), &
252 22 : ff_type%amb_info%bend_b(i), &
253 22 : ff_type%amb_info%bend_c(i), &
254 45 : k, theta0
255 : END DO
256 : ! Torsions
257 1 : WRITE (iw, 1003)
258 40 : DO i = 1, SIZE(ff_type%amb_info%torsion_a)
259 39 : k = cp_unit_from_cp2k(ff_type%amb_info%torsion_k(i), "kcalmol")
260 39 : m = ff_type%amb_info%torsion_m(i)
261 39 : phi0 = cp_unit_from_cp2k(ff_type%amb_info%torsion_phi0(i), "deg")
262 39 : WRITE (iw, 2003) ff_type%amb_info%torsion_a(i), &
263 39 : ff_type%amb_info%torsion_b(i), &
264 39 : ff_type%amb_info%torsion_c(i), &
265 39 : ff_type%amb_info%torsion_d(i), &
266 79 : k, m, phi0
267 : END DO
268 : ! Lennard-Jones
269 1 : WRITE (iw, 1005)
270 13 : DO i = 1, SIZE(ff_type%amb_info%nonbond_a)
271 12 : eps = cp_unit_from_cp2k(ff_type%amb_info%nonbond_eps(i), "kcalmol")
272 12 : sigma = cp_unit_from_cp2k(ff_type%amb_info%nonbond_rmin2(i), "angstrom")
273 12 : WRITE (iw, 2005) ff_type%amb_info%nonbond_a(i), &
274 25 : eps, sigma
275 : END DO
276 : END IF
277 : CASE (do_ff_g87, do_ff_g96)
278 0 : CPWARN("Dumping FF parameter file for GROMOS FF not implemented!")
279 : CASE (do_ff_undef)
280 2 : CPWARN("Dumping FF parameter file for INPUT FF not implemented!")
281 : END SELECT
282 2 : IF (iw > 0) THEN
283 1 : WRITE (iw, '(/,A)') "END"
284 : END IF
285 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
286 2 : "PRINT%FF_PARAMETER_FILE")
287 : END IF
288 2643 : CALL timestop(handle)
289 2643 : 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
|
