Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates QM/MM energy and forces
10 : !> \par History
11 : !> 2015 Factored out of force_env_methods.F
12 : !> \author Ole Schuett
13 : ! **************************************************************************************************
14 : MODULE qmmm_force
15 : USE cell_types, ONLY: cell_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_get_default_io_unit,&
18 : cp_logger_type
19 : USE cp_output_handling, ONLY: cp_iter_string,&
20 : cp_p_file,&
21 : cp_print_key_finished_output,&
22 : cp_print_key_should_output,&
23 : cp_print_key_unit_nr
24 : USE cp_result_methods, ONLY: cp_results_erase,&
25 : get_results,&
26 : put_results
27 : USE cp_result_types, ONLY: cp_result_type
28 : USE cp_subsys_types, ONLY: cp_subsys_get,&
29 : cp_subsys_type
30 : USE fist_energy_types, ONLY: fist_energy_type
31 : USE fist_environment_types, ONLY: fist_env_get
32 : USE fist_force, ONLY: fist_calc_energy_force
33 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
34 : section_vals_type
35 : USE kinds, ONLY: default_string_length,&
36 : dp
37 : USE particle_types, ONLY: particle_type
38 : USE physcon, ONLY: debye
39 : USE qmmm_gpw_energy, ONLY: qmmm_el_coupling
40 : USE qmmm_gpw_forces, ONLY: qmmm_forces
41 : USE qmmm_links_methods, ONLY: qmmm_added_chrg_coord,&
42 : qmmm_added_chrg_forces,&
43 : qmmm_link_Imomm_coord,&
44 : qmmm_link_Imomm_forces
45 : USE qmmm_types, ONLY: qmmm_env_type
46 : USE qmmm_types_low, ONLY: qmmm_links_type
47 : USE qmmm_util, ONLY: apply_qmmm_translate,&
48 : apply_qmmm_walls
49 : USE qs_energy_types, ONLY: qs_energy_type
50 : USE qs_environment_types, ONLY: get_qs_env
51 : USE qs_force, ONLY: qs_calc_energy_force
52 : USE qs_ks_qmmm_methods, ONLY: ks_qmmm_env_rebuild
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_force'
60 :
61 : PUBLIC :: qmmm_calc_energy_force
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief calculates the qm/mm energy and forces
67 : !> \param qmmm_env ...
68 : !> \param calc_force if also the forces should be calculated
69 : !> \param consistent_energies ...
70 : !> \param linres ...
71 : !> \par History
72 : !> 05.2004 created [fawzi]
73 : !> \author Fawzi Mohamed
74 : ! **************************************************************************************************
75 3802 : SUBROUTINE qmmm_calc_energy_force(qmmm_env, calc_force, consistent_energies, linres)
76 : TYPE(qmmm_env_type), POINTER :: qmmm_env
77 : LOGICAL, INTENT(IN) :: calc_force, consistent_energies, linres
78 :
79 : CHARACTER(LEN=default_string_length) :: description, iter
80 : INTEGER :: ip, j, nres, output_unit
81 3802 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index
82 : LOGICAL :: check, qmmm_added_chrg, qmmm_link, &
83 : qmmm_link_imomm
84 : REAL(KIND=dp) :: energy_mm, energy_qm
85 : REAL(KIND=dp), DIMENSION(3) :: dip_mm, dip_qm, dip_qmmm, max_coord, &
86 : min_coord
87 : TYPE(cell_type), POINTER :: mm_cell, qm_cell
88 : TYPE(cp_logger_type), POINTER :: logger
89 : TYPE(cp_result_type), POINTER :: results_mm, results_qm, results_qmmm
90 : TYPE(cp_subsys_type), POINTER :: subsys, subsys_mm, subsys_qm
91 : TYPE(fist_energy_type), POINTER :: fist_energy
92 3802 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm, particles_qm
93 : TYPE(qmmm_links_type), POINTER :: qmmm_links
94 : TYPE(qs_energy_type), POINTER :: qs_energy
95 : TYPE(section_vals_type), POINTER :: force_env_section, print_key
96 :
97 15208 : min_coord = HUGE(0.0_dp)
98 15208 : max_coord = -HUGE(0.0_dp)
99 3802 : qmmm_link = .FALSE.
100 3802 : qmmm_link_imomm = .FALSE.
101 3802 : qmmm_added_chrg = .FALSE.
102 3802 : logger => cp_get_default_logger()
103 3802 : output_unit = cp_logger_get_default_io_unit(logger)
104 3802 : NULLIFY (subsys_mm, subsys_qm, subsys, qm_atom_index, particles_mm, particles_qm, qm_cell, mm_cell)
105 3802 : NULLIFY (force_env_section, print_key, results_qmmm, results_qm, results_mm)
106 :
107 3802 : CALL get_qs_env(qmmm_env%qs_env, input=force_env_section)
108 3802 : print_key => section_vals_get_subs_vals(force_env_section, "QMMM%PRINT%DIPOLE")
109 :
110 3802 : CPASSERT(ASSOCIATED(qmmm_env))
111 :
112 3802 : CALL apply_qmmm_translate(qmmm_env)
113 :
114 3802 : CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
115 3802 : CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
116 3802 : qm_atom_index => qmmm_env%qm%qm_atom_index
117 3802 : qmmm_link = qmmm_env%qm%qmmm_link
118 3802 : qmmm_links => qmmm_env%qm%qmmm_links
119 3802 : qmmm_added_chrg = (qmmm_env%qm%move_mm_charges .OR. qmmm_env%qm%add_mm_charges)
120 3802 : IF (qmmm_link) THEN
121 280 : CPASSERT(ASSOCIATED(qmmm_links))
122 280 : IF (ASSOCIATED(qmmm_links%imomm)) qmmm_link_imomm = (SIZE(qmmm_links%imomm) /= 0)
123 : END IF
124 3802 : CPASSERT(ASSOCIATED(qm_atom_index))
125 :
126 3802 : particles_mm => subsys_mm%particles%els
127 3802 : particles_qm => subsys_qm%particles%els
128 :
129 15208 : DO j = 1, 3
130 11406 : IF (qm_cell%perd(j) == 1) CYCLE
131 63718 : DO ip = 1, SIZE(particles_qm)
132 : check = (DOT_PRODUCT(qm_cell%h_inv(j, :), particles_qm(ip)%r) >= 0.0) .AND. &
133 388944 : (DOT_PRODUCT(qm_cell%h_inv(j, :), particles_qm(ip)%r) <= 1.0)
134 48618 : IF (output_unit > 0 .AND. .NOT. check) THEN
135 0 : WRITE (unit=output_unit, fmt='("ip, j, pos, lat_pos ",2I6,6F12.5)') ip, j, &
136 0 : particles_qm(ip)%r, DOT_PRODUCT(qm_cell%h_inv(j, :), particles_qm(ip)%r)
137 : END IF
138 60024 : IF (.NOT. check) THEN
139 : CALL cp_abort(__LOCATION__, &
140 : "QM/MM QM atoms must be fully contained in the same image of the QM box "// &
141 0 : "- No wrapping of coordinates is allowed! ")
142 : END IF
143 : END DO
144 : END DO
145 :
146 : ! If present QM/MM links (just IMOMM) correct the position of the qm-link atom
147 3802 : IF (qmmm_link_imomm) CALL qmmm_link_Imomm_coord(qmmm_links, particles_qm, qm_atom_index)
148 :
149 : ! If add charges get their position NOW!
150 3802 : IF (qmmm_added_chrg) CALL qmmm_added_chrg_coord(qmmm_env%qm, particles_mm)
151 :
152 : ! Initialize ks_qmmm_env
153 3802 : CALL ks_qmmm_env_rebuild(qs_env=qmmm_env%qs_env, qmmm_env=qmmm_env%qm)
154 :
155 : ! Compute the short range QM/MM Electrostatic Potential
156 : CALL qmmm_el_coupling(qs_env=qmmm_env%qs_env, &
157 : qmmm_env=qmmm_env%qm, &
158 : mm_particles=particles_mm, &
159 3802 : mm_cell=mm_cell)
160 :
161 : ! Fist
162 3802 : CALL fist_calc_energy_force(qmmm_env%fist_env)
163 :
164 : ! Print Out information on fist energy calculation...
165 3802 : CALL fist_env_get(qmmm_env%fist_env, thermo=fist_energy)
166 3802 : energy_mm = fist_energy%pot
167 3802 : CALL cp_subsys_get(subsys_mm, results=results_mm)
168 :
169 : ! QS
170 3802 : CALL qs_calc_energy_force(qmmm_env%qs_env, calc_force, consistent_energies, linres)
171 :
172 : ! QM/MM Interaction Potential forces
173 : CALL qmmm_forces(qmmm_env%qs_env, &
174 : qmmm_env%qm, particles_mm, &
175 : mm_cell=mm_cell, &
176 3802 : calc_force=calc_force)
177 :
178 : ! Forces of quadratic wall on QM atoms
179 3802 : CALL apply_qmmm_walls(qmmm_env)
180 :
181 : ! Print Out information on QS energy calculation...
182 3802 : CALL get_qs_env(qmmm_env%qs_env, energy=qs_energy)
183 3802 : energy_qm = qs_energy%total
184 3802 : CALL cp_subsys_get(subsys_qm, results=results_qm)
185 :
186 : !TODO: is really results_qm == results_qmmm ???
187 3802 : CALL cp_subsys_get(subsys_qm, results=results_qmmm)
188 :
189 3802 : IF (calc_force) THEN
190 : ! If present QM/MM links (just IMOMM) correct the position of the qm-link atom
191 1744 : IF (qmmm_link_imomm) CALL qmmm_link_Imomm_forces(qmmm_links, particles_qm, qm_atom_index)
192 1744 : particles_mm => subsys_mm%particles%els
193 11288 : DO ip = 1, SIZE(qm_atom_index)
194 78096 : particles_mm(qm_atom_index(ip))%f = particles_mm(qm_atom_index(ip))%f + particles_qm(ip)%f
195 : END DO
196 : ! If add charges get rid of their derivatives right NOW!
197 1744 : IF (qmmm_added_chrg) CALL qmmm_added_chrg_forces(qmmm_env%qm, particles_mm)
198 : END IF
199 :
200 : ! Handle some output
201 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "QMMM%PRINT%DERIVATIVES", &
202 3802 : extension=".Log")
203 3802 : IF (output_unit > 0) THEN
204 0 : WRITE (unit=output_unit, fmt='(/1X,A,F15.9)') "Energy after QMMM calculation: ", energy_qm
205 0 : IF (calc_force) THEN
206 0 : WRITE (unit=output_unit, fmt='(/1X,A)') "Derivatives on all atoms after QMMM calculation: "
207 0 : DO ip = 1, SIZE(particles_mm)
208 0 : WRITE (unit=output_unit, fmt='(1X,3F15.9)') particles_mm(ip)%f
209 : END DO
210 : END IF
211 : END IF
212 : CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
213 3802 : "QMMM%PRINT%DERIVATIVES")
214 :
215 : ! Dipole
216 3802 : print_key => section_vals_get_subs_vals(force_env_section, "QMMM%PRINT%DIPOLE")
217 3802 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
218 : cp_p_file)) THEN
219 0 : description = '[DIPOLE]'
220 0 : CALL get_results(results=results_qm, description=description, n_rep=nres)
221 0 : CPASSERT(nres <= 1)
222 0 : CALL get_results(results=results_mm, description=description, n_rep=nres)
223 0 : CPASSERT(nres <= 1)
224 0 : CALL get_results(results=results_qm, description=description, values=dip_qm)
225 0 : CALL get_results(results=results_mm, description=description, values=dip_mm)
226 0 : dip_qmmm = dip_qm + dip_mm
227 0 : CALL cp_results_erase(results=results_qmmm, description=description)
228 0 : CALL put_results(results=results_qmmm, description=description, values=dip_qmmm)
229 :
230 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "QMMM%PRINT%DIPOLE", &
231 0 : extension=".Dipole")
232 0 : IF (output_unit > 0) THEN
233 0 : WRITE (unit=output_unit, fmt="(A)") "QMMM TOTAL DIPOLE"
234 : WRITE (unit=output_unit, fmt="(A,T31,A,T88,A)") &
235 0 : "# iter_level", "dipole(x,y,z)[atomic units]", &
236 0 : "dipole(x,y,z)[debye]"
237 0 : iter = cp_iter_string(logger%iter_info)
238 : WRITE (unit=output_unit, fmt="(a,6(es18.8))") &
239 0 : iter(1:15), dip_qmmm, dip_qmmm*debye
240 : END IF
241 : END IF
242 :
243 3802 : END SUBROUTINE qmmm_calc_energy_force
244 :
245 : END MODULE qmmm_force
|
