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