|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00011 MODULE qmmm_util 00012 USE cell_types, ONLY: cell_copy,& 00013 cell_type 00014 USE cp_subsys_types, ONLY: cp_subsys_type 00015 USE f77_blas 00016 USE force_env_types, ONLY: force_env_get,& 00017 force_env_type,& 00018 use_qmmm 00019 USE input_constants, ONLY: do_qmmm_wall_none,& 00020 do_qmmm_wall_quadratic,& 00021 do_qmmm_wall_reflective 00022 USE input_section_types, ONLY: section_vals_get,& 00023 section_vals_get_subs_vals,& 00024 section_vals_type,& 00025 section_vals_val_get 00026 USE kinds, ONLY: dp 00027 USE particle_types, ONLY: particle_type,& 00028 write_fist_particle_coordinates,& 00029 write_qs_particle_coordinates 00030 USE qmmm_types, ONLY: fist_subsys,& 00031 force_mixing_core_subsys,& 00032 force_mixing_extended_subsys,& 00033 primary_subsys,& 00034 qs_subsys 00035 USE qs_energy_types, ONLY: qs_energy_type 00036 USE qs_environment_types, ONLY: get_qs_env 00037 #include "cp_common_uses.h" 00038 00039 IMPLICIT NONE 00040 PRIVATE 00041 00042 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. 00043 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_util' 00044 PUBLIC :: apply_qmmm_walls_reflective,& 00045 apply_qmmm_walls,& 00046 apply_qmmm_translate,& 00047 apply_qmmm_translate_low,& 00048 spherical_cutoff_factor,& 00049 qmmm_force_mixing_active 00050 00051 CONTAINS 00052 00053 ! ***************************************************************************** 00062 SUBROUTINE apply_qmmm_walls(force_env, error) 00063 TYPE(force_env_type), POINTER :: force_env 00064 TYPE(cp_error_type), INTENT(inout) :: error 00065 00066 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls', 00067 routineP = moduleN//':'//routineN 00068 00069 INTEGER :: iwall_type 00070 LOGICAL :: explicit, failure 00071 TYPE(section_vals_type), POINTER :: walls_section 00072 00073 failure = .FALSE. 00074 walls_section => section_vals_get_subs_vals(force_env%root_section,& 00075 "FORCE_EVAL%QMMM%WALLS",error=error) 00076 CALL section_vals_get(walls_section, explicit=explicit, error=error) 00077 IF (explicit) THEN 00078 CALL section_vals_val_get(walls_section,"TYPE",i_val=iwall_type,error=error) 00079 SELECT CASE(iwall_type) 00080 CASE(do_qmmm_wall_quadratic) 00081 CALL apply_qmmm_walls_quadratic(force_env, walls_section, error) 00082 CASE(do_qmmm_wall_reflective) 00083 ! Do nothing.. reflective walls are applied directly in the integrator 00084 END SELECT 00085 ENDIF 00086 00087 END SUBROUTINE apply_qmmm_walls 00088 00089 ! ***************************************************************************** 00098 SUBROUTINE apply_qmmm_walls_reflective(force_env, error) 00099 TYPE(force_env_type), POINTER :: force_env 00100 TYPE(cp_error_type), INTENT(inout) :: error 00101 00102 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls_reflective', 00103 routineP = moduleN//':'//routineN 00104 00105 INTEGER :: ip, iwall_type, qm_index 00106 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 00107 LOGICAL :: explicit, failure, is_x(2), 00108 is_y(2), is_z(2) 00109 REAL(KIND=dp), DIMENSION(3) :: coord, qm_cell_diag, skin 00110 REAL(KIND=dp), DIMENSION(:), POINTER :: list 00111 TYPE(cell_type), POINTER :: mm_cell, qm_cell 00112 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 00113 TYPE(particle_type), DIMENSION(:), 00114 POINTER :: particles_mm 00115 TYPE(section_vals_type), POINTER :: walls_section 00116 00117 failure = .FALSE. 00118 SELECT CASE(force_env%in_use) 00119 CASE(use_qmmm) 00120 NULLIFY(subsys_mm, subsys_qm, qm_atom_index,particles_mm, qm_cell, mm_cell,& 00121 walls_section) 00122 walls_section => section_vals_get_subs_vals(force_env%root_section,"FORCE_EVAL%QMMM%WALLS",error=error) 00123 CALL section_vals_get(walls_section, explicit=explicit, error=error) 00124 IF (explicit) THEN 00125 NULLIFY(list) 00126 CALL section_vals_val_get(walls_section,"WALL_SKIN",r_vals=list,error=error) 00127 CALL section_vals_val_get(walls_section,"TYPE",i_val=iwall_type,error=error) 00128 skin(:) = list(:) 00129 ELSE 00130 ![NB] 00131 iwall_type=do_qmmm_wall_reflective 00132 skin(:) = 0.0_dp 00133 END IF 00134 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00135 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00136 CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure) 00137 CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure) 00138 00139 CALL cp_assert(iwall_type == do_qmmm_wall_none .OR. SIZE(force_env%sub_force_env) == 1, & 00140 cp_warning_level,cp_assertion_failed,routineP,& 00141 "Reflective walls for QM/MM are not implemented (or useful), when "//& 00142 "top level force_env has two or more sub_force_env allocated, which only happens when "//& 00143 "force mixing is active. Skipping!"//CPSourceFileRef) 00144 IF (SIZE(force_env%sub_force_env) /= 1) RETURN 00145 00146 CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& 00147 cell=mm_cell,subsys=subsys_mm,error=error) 00148 CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& 00149 cell=qm_cell,subsys=subsys_qm,error=error) 00150 qm_atom_index => force_env%qmmm_env%qm_atom_index 00151 CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) 00152 00153 qm_cell_diag = (/qm_cell%hmat(1,1),& 00154 qm_cell%hmat(2,2),& 00155 qm_cell%hmat(3,3)/) 00156 particles_mm => subsys_mm%particles%els 00157 DO ip=1,SIZE(qm_atom_index) 00158 qm_index = qm_atom_index(ip) 00159 coord = particles_mm(qm_index)%r 00160 IF (ANY(coord<skin).OR.ANY(coord>(qm_cell_diag-skin))) THEN 00161 IF (explicit) THEN 00162 IF (iwall_type==do_qmmm_wall_reflective) THEN 00163 ! Apply Walls 00164 is_x(1) = (coord(1)<skin(1)) 00165 is_x(2) = (coord(1)>(qm_cell_diag(1)-skin(1))) 00166 is_y(1) = (coord(2)<skin(2)) 00167 is_y(2) = (coord(2)>(qm_cell_diag(2)-skin(2))) 00168 is_z(1) = (coord(3)<skin(3)) 00169 is_z(2) = (coord(3)>(qm_cell_diag(3)-skin(3))) 00170 IF (ANY(is_x)) THEN 00171 ! X coordinate 00172 IF (is_x(1)) THEN 00173 particles_mm(qm_index)%v(1) = ABS(particles_mm(qm_index)%v(1)) 00174 ELSE IF (is_x(2)) THEN 00175 particles_mm(qm_index)%v(1) = -ABS(particles_mm(qm_index)%v(1)) 00176 END IF 00177 END IF 00178 IF (ANY(is_y)) THEN 00179 ! Y coordinate 00180 IF (is_y(1)) THEN 00181 particles_mm(qm_index)%v(2) = ABS(particles_mm(qm_index)%v(2)) 00182 ELSE IF (is_y(2)) THEN 00183 particles_mm(qm_index)%v(2) = -ABS(particles_mm(qm_index)%v(2)) 00184 END IF 00185 END IF 00186 IF (ANY(is_z)) THEN 00187 ! Z coordinate 00188 IF (is_z(1)) THEN 00189 particles_mm(qm_index)%v(3) = ABS(particles_mm(qm_index)%v(3)) 00190 ELSE IF (is_z(2)) THEN 00191 particles_mm(qm_index)%v(3) = -ABS(particles_mm(qm_index)%v(3)) 00192 END IF 00193 END IF 00194 ENDIF 00195 ELSE 00196 ! Otherwise print a warning and continue crossing cp2k's finger.. 00197 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00198 "One or few QM atoms are within the SKIN of the quantum box. Check your run "//& 00199 "and you may possibly consider: the activation of the QMMM WALLS "//& 00200 "around the QM box, switching ON the centering of the QM box or increase "//& 00201 "the size of the QM cell. CP2K CONTINUE but results could be meaningless. "//& 00202 CPSourceFileRef,& 00203 only_ionode=.TRUE.) 00204 END IF 00205 END IF 00206 END DO 00207 END SELECT 00208 00209 END SUBROUTINE apply_qmmm_walls_reflective 00210 00211 ! ***************************************************************************** 00220 SUBROUTINE apply_qmmm_walls_quadratic(force_env, walls_section, error) 00221 TYPE(force_env_type), POINTER :: force_env 00222 TYPE(section_vals_type), POINTER :: walls_section 00223 TYPE(cp_error_type), INTENT(inout) :: error 00224 00225 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls_quadratic', 00226 routineP = moduleN//':'//routineN 00227 00228 INTEGER :: ip, qm_index 00229 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 00230 LOGICAL :: failure, is_x(2), is_y(2), 00231 is_z(2) 00232 REAL(KIND=dp) :: k, wallenergy, wallforce 00233 REAL(KIND=dp), DIMENSION(3) :: coord, qm_cell_diag, skin 00234 REAL(KIND=dp), DIMENSION(:), POINTER :: list 00235 TYPE(cell_type), POINTER :: mm_cell, qm_cell 00236 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 00237 TYPE(particle_type), DIMENSION(:), 00238 POINTER :: particles_mm 00239 TYPE(qs_energy_type), POINTER :: energy 00240 00241 NULLIFY(list) 00242 CALL section_vals_val_get(walls_section,"WALL_SKIN",r_vals=list,error=error) 00243 CALL section_vals_val_get(walls_section,"K",r_val=k,error=error) 00244 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00245 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00246 CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure) 00247 CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure) 00248 00249 CALL cp_assert(SIZE(force_env%sub_force_env) == 1, & 00250 cp_warning_level,cp_assertion_failed,routineP,& 00251 "Quadratic walls for QM/MM are not implemented (or useful), when "//& 00252 "top level force_env has two or more sub_force_env allocated, which only happens when "//& 00253 "force mixing is active. Skipping!"//CPSourceFileRef) 00254 IF (SIZE(force_env%sub_force_env) /= 1) RETURN 00255 00256 CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& 00257 cell=mm_cell,subsys=subsys_mm,error=error) 00258 CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& 00259 cell=qm_cell,subsys=subsys_qm,error=error) 00260 00261 qm_atom_index => force_env%qmmm_env%qm_atom_index 00262 CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) 00263 00264 skin(:) = list(:) 00265 00266 qm_cell_diag = (/qm_cell%hmat(1,1),& 00267 qm_cell%hmat(2,2),& 00268 qm_cell%hmat(3,3)/) 00269 particles_mm => subsys_mm%particles%els 00270 wallenergy=0.0_dp 00271 DO ip=1,SIZE(qm_atom_index) 00272 qm_index = qm_atom_index(ip) 00273 coord = particles_mm(qm_index)%r 00274 IF (ANY(coord<skin).OR.ANY(coord>(qm_cell_diag-skin))) THEN 00275 is_x(1) = (coord(1)<skin(1)) 00276 is_x(2) = (coord(1)>(qm_cell_diag(1)-skin(1))) 00277 is_y(1) = (coord(2)<skin(2)) 00278 is_y(2) = (coord(2)>(qm_cell_diag(2)-skin(2))) 00279 is_z(1) = (coord(3)<skin(3)) 00280 is_z(2) = (coord(3)>(qm_cell_diag(3)-skin(3))) 00281 IF (is_x(1)) THEN 00282 wallforce=2.0_dp*k*(skin(1)-coord(1)) 00283 particles_mm(qm_index)%f(1)=particles_mm(qm_index)%f(1)+& 00284 wallforce 00285 wallenergy=wallenergy+wallforce*(skin(1)-coord(1))*0.5_dp 00286 ENDIF 00287 IF (is_x(2)) THEN 00288 wallforce=2.0_dp*k*((qm_cell_diag(1)-skin(1))-coord(1)) 00289 particles_mm(qm_index)%f(1)=particles_mm(qm_index)%f(1)+& 00290 wallforce 00291 wallenergy=wallenergy+wallforce*((qm_cell_diag(1)-skin(1))-& 00292 coord(1))*0.5_dp 00293 ENDIF 00294 IF (is_y(1)) THEN 00295 wallforce=2.0_dp*k*(skin(2)-coord(2)) 00296 particles_mm(qm_index)%f(2)=particles_mm(qm_index)%f(2)+& 00297 wallforce 00298 wallenergy=wallenergy+wallforce*(skin(2)-coord(2))*0.5_dp 00299 ENDIF 00300 IF (is_y(2)) THEN 00301 wallforce=2.0_dp*k*((qm_cell_diag(2)-skin(2))-coord(2)) 00302 particles_mm(qm_index)%f(2)=particles_mm(qm_index)%f(2)+& 00303 wallforce 00304 wallenergy=wallenergy+wallforce*((qm_cell_diag(2)-skin(2))-& 00305 coord(2))*0.5_dp 00306 ENDIF 00307 IF (is_z(1)) THEN 00308 wallforce=2.0_dp*k*(skin(3)-coord(3)) 00309 particles_mm(qm_index)%f(3)=particles_mm(qm_index)%f(3)+& 00310 wallforce 00311 wallenergy=wallenergy+wallforce*(skin(3)-coord(3))*0.5_dp 00312 ENDIF 00313 IF (is_z(2)) THEN 00314 wallforce=2.0_dp*k*((qm_cell_diag(3)-skin(3))-coord(3)) 00315 particles_mm(qm_index)%f(3)=particles_mm(qm_index)%f(3)+& 00316 wallforce 00317 wallenergy=wallenergy+wallforce*((qm_cell_diag(3)-skin(3))-& 00318 coord(3))*0.5_dp 00319 ENDIF 00320 ENDIF 00321 ENDDO 00322 00323 CALL get_qs_env(qs_env=force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%& 00324 force_env%qs_env, energy=energy,error=error) 00325 energy%total = energy%total + wallenergy 00326 00327 END SUBROUTINE apply_qmmm_walls_quadratic 00328 00329 ! ***************************************************************************** 00338 SUBROUTINE apply_qmmm_translate(force_env,error) 00339 TYPE(force_env_type), POINTER :: force_env 00340 TYPE(cp_error_type), INTENT(inout) :: error 00341 00342 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate', 00343 routineP = moduleN//':'//routineN 00344 00345 INTEGER :: ip 00346 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 00347 LOGICAL :: failure 00348 REAL(dp), POINTER :: charges(:) 00349 TYPE(cell_type), POINTER :: cell_core, cell_extended 00350 TYPE(cp_subsys_type), POINTER :: subsys_core, subsys_extended, 00351 subsys_fist, subsys_qs 00352 TYPE(particle_type), DIMENSION(:), 00353 POINTER :: particles_core, 00354 particles_extended, 00355 particles_mm, particles_qm 00356 TYPE(section_vals_type), POINTER :: subsys_section 00357 00358 SELECT CASE(force_env%in_use) 00359 00360 CASE(use_qmmm) 00361 IF (SIZE(force_env%sub_force_env) == 1) THEN 00362 CALL apply_qmmm_translate_low(force_env%sub_force_env(1)%force_env,error) 00363 ELSE IF (SIZE(force_env%sub_force_env) == 2) THEN 00364 ! want to center extended, and make core consistent with that 00365 CALL apply_qmmm_translate_low(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,error) 00366 00367 ! translate core fist particles 00368 CALL force_env_get(force_env%sub_force_env(force_mixing_extended_subsys)%force_env,& 00369 subsys=subsys_extended,cell=cell_extended,error=error) 00370 CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env,& 00371 subsys=subsys_core,cell=cell_core,error=error) 00372 particles_extended => subsys_extended%particles%els 00373 particles_core => subsys_core%particles%els 00374 DO ip=1,SIZE(particles_extended) 00375 particles_core(ip)%r = particles_extended(ip)%r 00376 END DO 00377 CALL cell_copy(cell_extended, cell_core, error) 00378 00379 ! also update core QM particles (like in apply_qmmm_translate_low) 00380 CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& 00381 subsys=subsys_fist,error=error) 00382 CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& 00383 subsys=subsys_qs,error=error) 00384 qm_atom_index => force_env%sub_force_env(force_mixing_core_subsys)%force_env%qmmm_env%qm_atom_index 00385 particles_qm => subsys_qs%particles%els 00386 particles_mm => subsys_fist%particles%els 00387 CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) 00388 DO ip=1,SIZE(qm_atom_index) 00389 particles_qm(ip)%r=particles_mm(qm_atom_index(ip))%r 00390 END DO 00391 00392 IF (debug_this_module) THEN 00393 CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& 00394 subsys=subsys_core,cell=cell_core,error=error) 00395 particles_core => subsys_core%particles%els 00396 subsys_section => section_vals_get_subs_vals(force_env%force_env_section, & 00397 "SUBSYS",error=error) 00398 CALL write_qs_particle_coordinates(particles_core,subsys_section,"QM/MM core calc first QM, then MM (0 charges)",error) 00399 CALL force_env_get(force_env%sub_force_env(force_mixing_core_subsys)%force_env%sub_force_env(fist_subsys)%force_env,& 00400 subsys=subsys_core,cell=cell_core,error=error) 00401 particles_core => subsys_core%particles%els 00402 ALLOCATE(charges(SIZE(particles_core))) 00403 charges = 0.0_dp 00404 CALL write_fist_particle_coordinates(particles_core,subsys_section,charges,error) 00405 DEALLOCATE(charges) 00406 ENDIF 00407 00408 ELSE 00409 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00410 ENDIF 00411 END SELECT 00412 END SUBROUTINE apply_qmmm_translate 00413 00414 ! ***************************************************************************** 00423 SUBROUTINE apply_qmmm_translate_low(force_env,error) 00424 TYPE(force_env_type), POINTER :: force_env 00425 TYPE(cp_error_type), INTENT(inout) :: error 00426 00427 CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate_low', 00428 routineP = moduleN//':'//routineN 00429 00430 INTEGER :: ip, unit_nr 00431 INTEGER, DIMENSION(:), POINTER :: qm_atom_index 00432 LOGICAL :: failure 00433 REAL(DP), POINTER :: charges(:) 00434 REAL(KIND=dp), DIMENSION(3) :: max_coord, min_coord, transl_v 00435 TYPE(cell_type), POINTER :: mm_cell, qm_cell 00436 TYPE(cp_logger_type), POINTER :: logger 00437 TYPE(cp_subsys_type), POINTER :: subsys_mm, subsys_qm 00438 TYPE(particle_type), DIMENSION(:), 00439 POINTER :: particles_mm, particles_qm 00440 TYPE(section_vals_type), POINTER :: subsys_section 00441 00442 SELECT CASE(force_env%in_use) 00443 CASE(use_qmmm) 00444 min_coord = HUGE(0.0_dp) 00445 max_coord = -HUGE(0.0_dp) 00446 failure = .FALSE. 00447 logger => cp_error_get_logger(error) 00448 NULLIFY(subsys_mm, subsys_qm, qm_atom_index, particles_mm, particles_qm,& 00449 subsys_section, qm_cell, mm_cell) 00450 00451 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00452 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00453 CPPrecondition(ASSOCIATED(force_env%qmmm_env),cp_failure_level,routineP,error,failure) 00454 CPPrecondition(force_env%qmmm_env%ref_count>0,cp_failure_level,routineP,error,failure) 00455 00456 CALL force_env_get(force_env%sub_force_env(fist_subsys)%force_env,& 00457 cell=mm_cell,subsys=subsys_mm,error=error) 00458 CALL force_env_get(force_env%sub_force_env(qs_subsys)%force_env,& 00459 cell=qm_cell,subsys=subsys_qm,error=error) 00460 qm_atom_index => force_env%qmmm_env%qm_atom_index 00461 CPPrecondition(ASSOCIATED(qm_atom_index),cp_failure_level,routineP,error,failure) 00462 00463 particles_qm => subsys_qm%particles%els 00464 particles_mm => subsys_mm%particles%els 00465 DO ip=1,SIZE(qm_atom_index) 00466 min_coord=MIN(min_coord,particles_mm(qm_atom_index(ip))%r) 00467 max_coord=MAX(max_coord,particles_mm(qm_atom_index(ip))%r) 00468 END DO 00469 IF (.NOT.force_env%qmmm_env%center_qm_subsys0) force_env%qmmm_env%do_translate = .FALSE. 00470 IF (force_env%qmmm_env%do_translate) THEN 00471 ! 00472 ! The first time we always translate all the system in order 00473 ! to centre the QM system in the box. 00474 ! 00475 transl_v = (max_coord + min_coord) / 2.0_dp 00476 transl_v(1) = transl_v(1) - qm_cell%hmat(1,1)/2.0_dp 00477 transl_v(2) = transl_v(2) - qm_cell%hmat(2,2)/2.0_dp 00478 transl_v(3) = transl_v(3) - qm_cell%hmat(3,3)/2.0_dp 00479 00480 IF (ANY(force_env%qmmm_env%utrasl /= 1.0_dp)) THEN 00481 transl_v = REAL( FLOOR(transl_v/force_env%qmmm_env%utrasl),KIND=dp) * 00482 force_env%qmmm_env%utrasl 00483 END IF 00484 force_env%qmmm_env%transl_v = force_env%qmmm_env%transl_v + transl_v 00485 particles_mm => subsys_mm%particles%els 00486 DO ip=1,subsys_mm%particles%n_els 00487 particles_mm(ip)%r = particles_mm(ip)%r - transl_v 00488 END DO 00489 unit_nr=cp_logger_get_default_io_unit(logger) 00490 IF (unit_nr>0) WRITE (unit=unit_nr,fmt='(/1X,A)')& 00491 " Translating the system in order to center the QM fragment in the QM box." 00492 IF (.NOT.force_env%qmmm_env%center_qm_subsys) force_env%qmmm_env%do_translate = .FALSE. 00493 END IF 00494 particles_mm => subsys_mm%particles%els 00495 DO ip=1,SIZE(qm_atom_index) 00496 particles_qm(ip)%r=particles_mm(qm_atom_index(ip))%r 00497 END DO 00498 00499 subsys_section => section_vals_get_subs_vals(force_env%force_env_section, & 00500 "SUBSYS",error=error) 00501 CALL write_qs_particle_coordinates(particles_qm,subsys_section,"QM/MM first QM, then MM (0 charges)",error) 00502 ALLOCATE(charges(SIZE(particles_mm))) 00503 charges = 0.0_dp 00504 CALL write_fist_particle_coordinates(particles_mm,subsys_section,charges,error) 00505 DEALLOCATE(charges) 00506 00507 END SELECT 00508 00509 END SUBROUTINE apply_qmmm_translate_low 00510 00511 ! ***************************************************************************** 00519 SUBROUTINE spherical_cutoff_factor(spherical_cutoff, rij, factor, error) 00520 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: spherical_cutoff 00521 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij 00522 REAL(KIND=dp), INTENT(OUT) :: factor 00523 TYPE(cp_error_type), INTENT(inout) :: error 00524 00525 CHARACTER(len=*), PARAMETER :: routineN = 'spherical_cutoff_factor', 00526 routineP = moduleN//':'//routineN 00527 00528 REAL(KIND=dp) :: r, r0 00529 00530 r = SQRT(DOT_PRODUCT(rij,rij)) 00531 r0 = spherical_cutoff(1)-20.0_dp*spherical_cutoff(2) 00532 factor = 0.5_dp*(1.0_dp-TANH((r-r0)/spherical_cutoff(2))) 00533 00534 END SUBROUTINE spherical_cutoff_factor 00535 00536 RECURSIVE FUNCTION qmmm_force_mixing_active(force_env, error) RESULT(active) 00537 TYPE(force_env_type), POINTER :: force_env 00538 TYPE(cp_error_type), INTENT(inout) :: error 00539 LOGICAL :: active 00540 00541 INTEGER :: iforce_eval, nforce_eval 00542 TYPE(section_vals_type), POINTER :: qmmm_force_mixing_section 00543 00544 active = .FALSE. 00545 qmmm_force_mixing_section => section_vals_get_subs_vals(force_env%force_env_section,"QMMM%FORCE_MIXING",& 00546 can_return_null=.TRUE.,error=error) 00547 IF (ASSOCIATED(qmmm_force_mixing_section)) CALL section_vals_get(qmmm_force_mixing_section,& 00548 explicit=active,error=error) 00549 00550 IF (active) RETURN 00551 00552 ! top level wasn't a QMMM with FORCE_MIXING, so loop over sub force envs, e.g. for mixed force env 00553 nforce_eval = SIZE(force_env%sub_force_env) 00554 DO iforce_eval=1, nforce_eval 00555 ! if the current force env is QMMM, then any sub force envs are dups of it generated by FORCE_MIXING, so explicitly avoid recursing into those 00556 IF (force_env%in_use /= use_qmmm) THEN 00557 active = qmmm_force_mixing_active(force_env%sub_force_env(iforce_eval)%force_env, error) 00558 IF (active) RETURN 00559 ENDIF 00560 END DO 00561 00562 END FUNCTION qmmm_force_mixing_active 00563 00564 END MODULE qmmm_util
1.7.3