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 : !> \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 48618 : IF (.NOT. check) &
139 : CALL cp_abort(__LOCATION__, &
140 : "QM/MM QM atoms must be fully contained in the same image of the QM box "// &
141 11406 : "- No wrapping of coordinates is allowed! ")
142 : END DO
143 : END DO
144 :
145 : ! If present QM/MM links (just IMOMM) correct the position of the qm-link atom
146 3802 : IF (qmmm_link_imomm) CALL qmmm_link_Imomm_coord(qmmm_links, particles_qm, qm_atom_index)
147 :
148 : ! If add charges get their position NOW!
149 3802 : IF (qmmm_added_chrg) CALL qmmm_added_chrg_coord(qmmm_env%qm, particles_mm)
150 :
151 : ! Initialize ks_qmmm_env
152 3802 : CALL ks_qmmm_env_rebuild(qs_env=qmmm_env%qs_env, qmmm_env=qmmm_env%qm)
153 :
154 : ! Compute the short range QM/MM Electrostatic Potential
155 : CALL qmmm_el_coupling(qs_env=qmmm_env%qs_env, &
156 : qmmm_env=qmmm_env%qm, &
157 : mm_particles=particles_mm, &
158 3802 : mm_cell=mm_cell)
159 :
160 : ! Fist
161 3802 : CALL fist_calc_energy_force(qmmm_env%fist_env)
162 :
163 : ! Print Out information on fist energy calculation...
164 3802 : CALL fist_env_get(qmmm_env%fist_env, thermo=fist_energy)
165 3802 : energy_mm = fist_energy%pot
166 3802 : CALL cp_subsys_get(subsys_mm, results=results_mm)
167 :
168 : ! QS
169 3802 : CALL qs_calc_energy_force(qmmm_env%qs_env, calc_force, consistent_energies, linres)
170 :
171 : ! QM/MM Interaction Potential forces
172 : CALL qmmm_forces(qmmm_env%qs_env, &
173 : qmmm_env%qm, particles_mm, &
174 : mm_cell=mm_cell, &
175 3802 : calc_force=calc_force)
176 :
177 : ! Forces of quadratic wall on QM atoms
178 3802 : CALL apply_qmmm_walls(qmmm_env)
179 :
180 : ! Print Out information on QS energy calculation...
181 3802 : CALL get_qs_env(qmmm_env%qs_env, energy=qs_energy)
182 3802 : energy_qm = qs_energy%total
183 3802 : CALL cp_subsys_get(subsys_qm, results=results_qm)
184 :
185 : !TODO: is really results_qm == results_qmmm ???
186 3802 : CALL cp_subsys_get(subsys_qm, results=results_qmmm)
187 :
188 3802 : IF (calc_force) THEN
189 : ! If present QM/MM links (just IMOMM) correct the position of the qm-link atom
190 1744 : IF (qmmm_link_imomm) CALL qmmm_link_Imomm_forces(qmmm_links, particles_qm, qm_atom_index)
191 1744 : particles_mm => subsys_mm%particles%els
192 11288 : DO ip = 1, SIZE(qm_atom_index)
193 78096 : particles_mm(qm_atom_index(ip))%f = particles_mm(qm_atom_index(ip))%f + particles_qm(ip)%f
194 : END DO
195 : ! If add charges get rid of their derivatives right NOW!
196 1744 : IF (qmmm_added_chrg) CALL qmmm_added_chrg_forces(qmmm_env%qm, particles_mm)
197 : END IF
198 :
199 : ! Handle some output
200 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "QMMM%PRINT%DERIVATIVES", &
201 3802 : extension=".Log")
202 3802 : IF (output_unit > 0) THEN
203 0 : WRITE (unit=output_unit, fmt='(/1X,A,F15.9)') "Energy after QMMM calculation: ", energy_qm
204 0 : IF (calc_force) THEN
205 0 : WRITE (unit=output_unit, fmt='(/1X,A)') "Derivatives on all atoms after QMMM calculation: "
206 0 : DO ip = 1, SIZE(particles_mm)
207 0 : WRITE (unit=output_unit, fmt='(1X,3F15.9)') particles_mm(ip)%f
208 : END DO
209 : END IF
210 : END IF
211 : CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
212 3802 : "QMMM%PRINT%DERIVATIVES")
213 :
214 : ! Dipole
215 3802 : print_key => section_vals_get_subs_vals(force_env_section, "QMMM%PRINT%DIPOLE")
216 3802 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
217 : cp_p_file)) THEN
218 0 : description = '[DIPOLE]'
219 0 : CALL get_results(results=results_qm, description=description, n_rep=nres)
220 0 : CPASSERT(nres <= 1)
221 0 : CALL get_results(results=results_mm, description=description, n_rep=nres)
222 0 : CPASSERT(nres <= 1)
223 0 : CALL get_results(results=results_qm, description=description, values=dip_qm)
224 0 : CALL get_results(results=results_mm, description=description, values=dip_mm)
225 0 : dip_qmmm = dip_qm + dip_mm
226 0 : CALL cp_results_erase(results=results_qmmm, description=description)
227 0 : CALL put_results(results=results_qmmm, description=description, values=dip_qmmm)
228 :
229 : output_unit = cp_print_key_unit_nr(logger, force_env_section, "QMMM%PRINT%DIPOLE", &
230 0 : extension=".Dipole")
231 0 : IF (output_unit > 0) THEN
232 0 : WRITE (unit=output_unit, fmt="(A)") "QMMM TOTAL DIPOLE"
233 : WRITE (unit=output_unit, fmt="(A,T31,A,T88,A)") &
234 0 : "# iter_level", "dipole(x,y,z)[atomic units]", &
235 0 : "dipole(x,y,z)[debye]"
236 0 : iter = cp_iter_string(logger%iter_info)
237 : WRITE (unit=output_unit, fmt="(a,6(es18.8))") &
238 0 : iter(1:15), dip_qmmm, dip_qmmm*debye
239 : END IF
240 : END IF
241 :
242 3802 : END SUBROUTINE qmmm_calc_energy_force
243 :
244 : END MODULE qmmm_force
|
