1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines to compute energy and forces in a QM/MM calculation 8!> \par History 9!> 05.2004 created [tlaino] 10!> \author Teodoro Laino 11! ************************************************************************************************** 12MODULE qmmm_gpw_forces 13 USE cell_types, ONLY: cell_type,& 14 pbc 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_log_handling, ONLY: cp_get_default_logger,& 17 cp_logger_type 18 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 19 cp_print_key_unit_nr 20 USE cp_para_types, ONLY: cp_para_env_type 21 USE cp_spline_utils, ONLY: pw_restrict_s3,& 22 spline3_nopbc_interp,& 23 spline3_pbc_interp 24 USE cube_utils, ONLY: cube_info_type 25 USE input_constants, ONLY: do_par_atom,& 26 do_qmmm_coulomb,& 27 do_qmmm_gauss,& 28 do_qmmm_none,& 29 do_qmmm_pcharge,& 30 do_qmmm_swave 31 USE input_section_types, ONLY: section_vals_get_subs_vals,& 32 section_vals_type,& 33 section_vals_val_get 34 USE kinds, ONLY: dp 35 USE message_passing, ONLY: mp_irecv,& 36 mp_isend,& 37 mp_sum,& 38 mp_wait 39 USE mm_collocate_potential, ONLY: collocate_gf_rspace_NoPBC,& 40 integrate_gf_rspace_NoPBC 41 USE particle_types, ONLY: particle_type 42 USE pw_env_types, ONLY: pw_env_get,& 43 pw_env_type 44 USE pw_methods, ONLY: pw_axpy,& 45 pw_integral_ab,& 46 pw_transfer,& 47 pw_zero 48 USE pw_pool_types, ONLY: pw_pool_create_pw,& 49 pw_pool_give_back_pw,& 50 pw_pool_p_type,& 51 pw_pool_type,& 52 pw_pools_create_pws,& 53 pw_pools_give_back_pws 54 USE pw_types, ONLY: REALDATA3D,& 55 REALSPACE,& 56 pw_p_type,& 57 pw_type 58 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,& 59 qmmm_gaussian_type 60 USE qmmm_gpw_energy, ONLY: qmmm_elec_with_gaussian,& 61 qmmm_elec_with_gaussian_LG,& 62 qmmm_elec_with_gaussian_LR 63 USE qmmm_se_forces, ONLY: deriv_se_qmmm_matrix 64 USE qmmm_types_low, ONLY: qmmm_env_qm_type,& 65 qmmm_per_pot_p_type,& 66 qmmm_per_pot_type,& 67 qmmm_pot_p_type,& 68 qmmm_pot_type 69 USE qmmm_util, ONLY: spherical_cutoff_factor 70 USE qmmm_tb_methods, ONLY: deriv_tb_qmmm_matrix,& 71 deriv_tb_qmmm_matrix_pc 72 USE qs_energy_types, ONLY: qs_energy_type 73 USE qs_environment_types, ONLY: get_qs_env,& 74 qs_environment_type 75 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type 76 USE qs_rho_types, ONLY: qs_rho_get,& 77 qs_rho_type 78#include "./base/base_uses.f90" 79 80 IMPLICIT NONE 81 82 PRIVATE 83 LOGICAL, PARAMETER, PRIVATE :: debug_this_module=.FALSE. 84 REAL(KIND=dp), PARAMETER, PRIVATE :: Dx = 0.01_dp ! Debug Variables 85 REAL(KIND=dp), PARAMETER, PRIVATE :: MaxErr = 10.0_dp ! Debug Variables 86 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_forces' 87 PUBLIC :: qmmm_forces 88 89CONTAINS 90 91! ************************************************************************************************** 92!> \brief General driver to Compute the contribution 93!> to the forces due to the QM/MM potential 94!> \param qs_env ... 95!> \param qmmm_env ... 96!> \param mm_particles ... 97!> \param calc_force ... 98!> \param mm_cell ... 99!> \par History 100!> 06.2004 created [tlaino] 101!> \author Teodoro Laino 102! ************************************************************************************************** 103 SUBROUTINE qmmm_forces(qs_env,qmmm_env,mm_particles,calc_force,mm_cell) 104 TYPE(qs_environment_type), POINTER :: qs_env 105 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 106 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 107 LOGICAL, INTENT(in), OPTIONAL :: calc_force 108 TYPE(cell_type), POINTER :: mm_cell 109 110 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces', routineP = moduleN//':'//routineN 111 112 INTEGER :: handle, iatom, image_IndMM, Imm, IndMM, & 113 ispin, iw 114 LOGICAL :: gapw, need_f, periodic 115 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges, & 116 Forces_added_shells 117 TYPE(cp_logger_type), POINTER :: logger 118 TYPE(cp_para_env_type), POINTER :: para_env 119 TYPE(dft_control_type), POINTER :: dft_control 120 TYPE(pw_env_type), POINTER :: pw_env 121 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 122 TYPE(pw_p_type), POINTER :: rho0_s_gs, rho_core 123 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 124 TYPE(pw_pool_type), POINTER :: auxbas_pool 125 TYPE(pw_type), POINTER :: rho_tot_r, rho_tot_r2 126 TYPE(qs_energy_type), POINTER :: energy 127 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc 128 TYPE(qs_rho_type), POINTER :: rho 129 TYPE(section_vals_type), POINTER :: input_section, interp_section, & 130 print_section 131 132 CALL timeset(routineN,handle) 133 need_f = .TRUE. 134 periodic = qmmm_env%periodic 135 IF (PRESENT(calc_force)) need_f = calc_force 136 NULLIFY(dft_control, ks_qmmm_env_loc, rho, pw_env, rho_tot_r, energy, Forces,& 137 Forces_added_charges,input_section,rho0_s_gs, rho_r) 138 CALL get_qs_env(qs_env=qs_env,& 139 rho=rho,& 140 rho_core=rho_core,& 141 pw_env=pw_env,& 142 energy=energy,& 143 para_env=para_env,& 144 input=input_section,& 145 rho0_s_gs=rho0_s_gs,& 146 dft_control=dft_control) 147 148 CALL qs_rho_get(rho, rho_r=rho_r) 149 150 logger => cp_get_default_logger() 151 ks_qmmm_env_loc => qs_env%ks_qmmm_env 152 interp_section => section_vals_get_subs_vals(input_section,"QMMM%INTERPOLATOR") 153 print_section => section_vals_get_subs_vals(input_section,"QMMM%PRINT") 154 iw=cp_print_key_unit_nr(logger,print_section,"PROGRAM_RUN_INFO",& 155 extension=".qmmmLog") 156 gapw = dft_control%qs_control%gapw 157 ! If forces are required allocate these temporary arrays 158 IF (need_f) THEN 159 ALLOCATE(Forces(3,qmmm_env%num_mm_atoms)) 160 ALLOCATE(Forces_added_charges(3,qmmm_env%added_charges%num_mm_atoms)) 161 ALLOCATE(Forces_added_shells(3,qmmm_env%added_shells%num_mm_atoms)) 162 Forces(:,:) = 0.0_dp 163 Forces_added_charges(:,:) = 0.0_dp 164 Forces_added_shells(:,:) = 0.0_dp 165 END IF 166 IF (dft_control%qs_control%semi_empirical) THEN 167 ! SEMIEMPIRICAL 168 SELECT CASE(qmmm_env%qmmm_coupl_type) 169 CASE(do_qmmm_coulomb) 170 CALL deriv_se_qmmm_matrix(qs_env,qmmm_env,mm_particles,mm_cell,para_env,& 171 need_f,Forces,Forces_added_charges) 172 CASE(do_qmmm_pcharge) 173 CPABORT("Point Charge QM/MM electrostatic coupling not yet implemented for SE.") 174 CASE(do_qmmm_gauss,do_qmmm_swave) 175 CPABORT("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.") 176 CASE(do_qmmm_none) 177 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')& 178 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!" 179 CASE DEFAULT 180 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')"Unknown Coupling..." 181 CPABORT("") 182 END SELECT 183 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN 184 ! DFTB 185 SELECT CASE(qmmm_env%qmmm_coupl_type) 186 CASE(do_qmmm_none) 187 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')& 188 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!" 189 CASE(do_qmmm_coulomb) 190 CALL deriv_tb_qmmm_matrix(qs_env,qmmm_env,mm_particles,mm_cell,para_env,& 191 need_f,Forces,Forces_added_charges) 192 CASE(do_qmmm_pcharge) 193 CALL deriv_tb_qmmm_matrix_pc(qs_env,qmmm_env,mm_particles,mm_cell,para_env,& 194 need_f,Forces,Forces_added_charges) 195 CASE(do_qmmm_gauss,do_qmmm_swave) 196 CPABORT("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for DFTB.") 197 CASE DEFAULT 198 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')"Unknown Coupling..." 199 CPABORT("") 200 END SELECT 201 IF (need_f) THEN 202 Forces(:,:) = Forces(:,:) / REAL(para_env%num_pe,KIND=dp) 203 Forces_added_charges(:,:) = Forces_added_charges(:,:) / REAL(para_env%num_pe,KIND=dp) 204 ENDIF 205 ELSE 206 ! GPW/GAPW 207 CALL pw_env_get(pw_env=pw_env,& 208 pw_pools=pw_pools,& 209 auxbas_pw_pool=auxbas_pool) 210 CALL pw_pool_create_pw(auxbas_pool,rho_tot_r,& 211 in_space=REALSPACE,& 212 use_data=REALDATA3D) 213 ! IF GAPW the core charge is replaced by the compensation charge 214 IF(gapw) THEN 215 IF( dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN 216 CALL pw_transfer(rho_core%pw,rho_tot_r) 217 energy%qmmm_nu = pw_integral_ab ( rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace%pw) 218 CALL pw_pool_create_pw(auxbas_pool,rho_tot_r2,& 219 in_space=REALSPACE,& 220 use_data=REALDATA3D) 221 CALL pw_transfer(rho0_s_gs%pw,rho_tot_r2) 222 CALL pw_axpy(rho_tot_r2,rho_tot_r) 223 CALL pw_pool_give_back_pw(auxbas_pool,rho_tot_r2,accept_non_compatible=.TRUE.) 224 ELSE 225 CALL pw_transfer(rho0_s_gs%pw,rho_tot_r) 226 ! 227 ! QM/MM Nuclear Electrostatic Potential already included through rho0 228 ! 229 energy%qmmm_nu = 0.0_dp 230 END IF 231 ELSE 232 CALL pw_transfer(rho_core%pw,rho_tot_r) 233 ! 234 ! Computes the QM/MM Nuclear Electrostatic Potential 235 ! 236 energy%qmmm_nu = pw_integral_ab ( rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace%pw) 237 END IF 238 IF (need_f) THEN 239 ! 240 DO ispin=1,SIZE(rho_r) 241 CALL pw_axpy(rho_r(ispin)%pw,rho_tot_r) 242 END DO 243 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')"Evaluating forces on MM atoms due to the:" 244 ! Electrostatic Interaction type... 245 SELECT CASE(qmmm_env%qmmm_coupl_type) 246 CASE(do_qmmm_coulomb) 247 CPABORT("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.") 248 CASE(do_qmmm_pcharge) 249 CPABORT("Point Charge QM/MM electrostatic coupling not yet implemented for GPW/GAPW.") 250 CASE(do_qmmm_gauss,do_qmmm_swave) 251 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')& 252 "- QM/MM Coupling computed collocating the Gaussian Potential Functions." 253 CALL qmmm_forces_with_gaussian(rho=rho_tot_r,& 254 qmmm_env=qmmm_env,& 255 mm_particles=mm_particles,& 256 aug_pools=qmmm_env%aug_pools,& 257 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid,& 258 coarser_grid=qmmm_env%gridlevel_info%coarser_grid,& 259 para_env=para_env,& 260 pw_pools=pw_pools,& 261 eps_mm_rspace=qmmm_env%eps_mm_rspace,& 262 cube_info=ks_qmmm_env_loc%cube_info,& 263 Forces=Forces,& 264 Forces_added_charges=Forces_added_charges,& 265 Forces_added_shells=Forces_added_shells,& 266 interp_section=interp_section,& 267 iw=iw,& 268 mm_cell=mm_cell) 269 CASE(do_qmmm_none) 270 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')& 271 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!" 272 CASE DEFAULT 273 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')"- Unknown Coupling..." 274 CPABORT("") 275 END SELECT 276 END IF 277 END IF 278 ! Correct Total Energy adding the contribution of the QM/MM nuclear interaction 279 energy%total = energy%total + energy%qmmm_nu 280 ! Proceed if gradients are requested.. 281 IF (need_f) THEN 282 !ikuo Temporary change to alleviate compiler problems on Intel with 283 !array dimension of 0 284 IF(qmmm_env%num_mm_atoms/=0) CALL mp_sum(Forces,para_env%group) 285 IF(qmmm_env%added_charges%num_mm_atoms/=0) CALL mp_sum(Forces_added_charges,para_env%group) 286 IF(qmmm_env%added_shells%num_mm_atoms/=0) CALL mp_sum(Forces_added_shells,para_env%group) 287 ! Debug Forces 288 IF (debug_this_module) THEN 289 IF ( dft_control%qs_control%semi_empirical.OR.& 290 dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)THEN 291 WRITE(iw,*)"NO DEBUG AVAILABLE in module"//TRIM(routineN) 292 ELSE 293 ! Print Out Forces 294 IF (iw>0) THEN 295 DO Imm = 1, SIZE(qmmm_env%mm_atom_index) 296 WRITE(iw,*)"ANALYTICAL FORCES:" 297 IndMM = qmmm_env%mm_atom_index(Imm) 298 WRITE(iw,'(I6,3F15.9)') IndMM, Forces(:,Imm) 299 END DO 300 END IF 301 CALL qmmm_debug_forces(rho=rho_tot_r,& 302 qs_env=qs_env,& 303 qmmm_env=qmmm_env,& 304 Analytical_Forces=Forces,& 305 mm_particles=mm_particles,& 306 mm_atom_index=qmmm_env%mm_atom_index,& 307 num_mm_atoms=qmmm_env%num_mm_atoms,& 308 interp_section=interp_section,& 309 mm_cell=mm_cell) 310 END IF 311 ENDIF 312 END IF 313 ! Give back rho_tot_t to auxbas_pool only for GPW/GAPW 314 IF ((.NOT.dft_control%qs_control%semi_empirical).AND.& 315 (.NOT.dft_control%qs_control%dftb).AND.(.NOT.dft_control%qs_control%xtb)) THEN 316 CALL pw_pool_give_back_pw(auxbas_pool,rho_tot_r, accept_non_compatible=.TRUE.) 317 END IF 318 IF (iw>0) THEN 319 IF(.NOT. gapw) WRITE(iw,'(T2,"QMMM|",1X,A,T66,F15.9)')& 320 "QM/MM Nuclear Electrostatic Potential :",energy%qmmm_nu 321 WRITE(iw,'(T2,"QMMM|",1X,A,T66,F15.9)')& 322 "QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):",energy%total 323 WRITE(iw,'(T2,"QMMM|",1X,A)')"MM energy NOT included in the above term!"//& 324 " Check for: FORCE_EVAL ( QMMM )" 325 WRITE(iw,'(T2,"QMMM|",1X,A)')"that includes both QM, QMMM and MM energy terms!" 326 END IF 327 IF (need_f) THEN 328 ! Transfer Forces 329 DO Imm = 1, qmmm_env%num_mm_atoms 330 IndMM = qmmm_env%mm_atom_index(Imm) 331 332 !add image forces to Forces 333 IF(qmmm_env%image_charge) THEN 334 DO iatom=1,qmmm_env%num_image_mm_atoms 335 image_IndMM=qmmm_env%image_charge_pot%image_mm_list(iatom) 336 IF(image_IndMM.eq.IndMM) THEN 337 Forces(:,Imm)=Forces(:,Imm)& 338 +qmmm_env%image_charge_pot%image_forcesMM(:,iatom) 339 ENDIF 340 ENDDO 341 ENDIF 342 343 ! Hack: In Forces there the gradients indeed... 344 ! Minux sign to take care of this misunderstanding... 345 mm_particles(IndMM)%f(:) = - Forces(:,Imm) + mm_particles(IndMM)%f(:) 346 END DO 347 DEALLOCATE(Forces) 348 IF (qmmm_env%move_mm_charges.OR.qmmm_env%add_mm_charges) THEN 349 DO Imm = 1, qmmm_env%added_charges%num_mm_atoms 350 IndMM = qmmm_env%added_charges%mm_atom_index(Imm) 351 ! Hack: In Forces there the gradients indeed... 352 ! Minux sign to take care of this misunderstanding... 353 qmmm_env%added_charges%added_particles(IndMM)%f(:) = - Forces_added_charges(:,Imm) 354 END DO 355 END IF 356 DEALLOCATE(Forces_added_charges) 357 IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN 358 DO Imm = 1, qmmm_env%added_shells%num_mm_atoms 359 IndMM = qmmm_env%added_shells%mm_core_index(Imm) 360 ! Hack: In Forces there the gradients indeed... 361 ! Minux sign to take care of this misunderstanding... 362 qmmm_env%added_shells%added_particles(Imm)%f(:) = qmmm_env%added_shells%added_particles(Imm)%f(:) - & 363 Forces_added_shells(:,Imm) 364 365 END DO 366 END IF 367 DEALLOCATE(Forces_added_shells) 368 END IF 369 CALL cp_print_key_finished_output(iw,logger,print_section,"PROGRAM_RUN_INFO") 370 CALL timestop(handle) 371 372 END SUBROUTINE qmmm_forces 373 374! ************************************************************************************************** 375!> \brief Evaluates the contribution to the forces due to the 376!> QM/MM potential computed collocating the Electrostatic 377!> Gaussian Potential. 378!> \param rho ... 379!> \param qmmm_env ... 380!> \param mm_particles ... 381!> \param aug_pools ... 382!> \param auxbas_grid ... 383!> \param coarser_grid ... 384!> \param cube_info ... 385!> \param para_env ... 386!> \param eps_mm_rspace ... 387!> \param pw_pools ... 388!> \param Forces ... 389!> \param Forces_added_charges ... 390!> \param Forces_added_shells ... 391!> \param interp_section ... 392!> \param iw ... 393!> \param mm_cell ... 394!> \par History 395!> 06.2004 created [tlaino] 396!> \author Teodoro Laino 397! ************************************************************************************************** 398 SUBROUTINE qmmm_forces_with_gaussian(rho, qmmm_env, mm_particles, & 399 aug_pools, auxbas_grid, coarser_grid, cube_info, para_env, & 400 eps_mm_rspace, pw_pools, Forces, Forces_added_charges, Forces_added_shells,& 401 interp_section, iw, mm_cell) 402 TYPE(pw_type), POINTER :: rho 403 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 404 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 405 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 406 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid 407 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info 408 TYPE(cp_para_env_type), POINTER :: para_env 409 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace 410 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 411 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges, & 412 Forces_added_shells 413 TYPE(section_vals_type), POINTER :: interp_section 414 INTEGER, INTENT(IN) :: iw 415 TYPE(cell_type), POINTER :: mm_cell 416 417 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces_with_gaussian', & 418 routineP = moduleN//':'//routineN 419 420 INTEGER :: group, handle, i, igrid, j, k, & 421 kind_interp, me, ngrids, request, stat 422 INTEGER, DIMENSION(3) :: glb, gub, lb, ub 423 INTEGER, DIMENSION(:), POINTER :: pos_of_x 424 LOGICAL :: shells 425 REAL(KIND=dp), DIMENSION(:, :), POINTER :: tmp 426 TYPE(pw_p_type), DIMENSION(:), POINTER :: grids 427 428! Statements 429 430 CALL timeset(routineN,handle) 431 NULLIFY(grids,tmp) 432 CPASSERT(ASSOCIATED(mm_particles)) 433 CPASSERT(ASSOCIATED(qmmm_env%mm_atom_chrg)) 434 CPASSERT(ASSOCIATED(qmmm_env%mm_atom_index)) 435 CPASSERT(ASSOCIATED(Forces)) 436 !Statements 437 ngrids=SIZE(pw_pools) 438 CALL pw_pools_create_pws(aug_pools,grids,use_data=REALDATA3D,in_space=REALSPACE) 439 DO igrid=1,ngrids 440 CALL pw_zero(grids(igrid)%pw) 441 END DO 442 ! Collocate Density on multigrids 443 lb = rho%pw_grid%bounds_local(1,:) 444 ub = rho%pw_grid%bounds_local(2,:) 445 grids(auxbas_grid) % pw % cr3d(lb(1):ub(1),& 446 lb(2):ub(2),& 447 lb(3):ub(3) )= rho % cr3d 448 ! copy the boundaries 449 DO i=lb(1),ub(1) 450 grids(auxbas_grid) % pw %cr3d(i,ub(2)+1,ub(3)+1)=rho%cr3d(i,lb(2),lb(3)) 451 END DO 452 DO k=lb(3),ub(3) 453 DO i=lb(1),ub(1) 454 grids(auxbas_grid) % pw %cr3d(i,ub(2)+1,k)=rho%cr3d(i,lb(2),k) 455 END DO 456 END DO 457 DO j=lb(2),ub(2) 458 DO i=lb(1),ub(1) 459 grids(auxbas_grid) % pw %cr3d(i,j,ub(3)+1)=rho%cr3d(i,j,lb(3)) 460 END DO 461 END DO 462 pos_of_x => grids(auxbas_grid)%pw%pw_grid%para%pos_of_x 463 group=grids(auxbas_grid)%pw%pw_grid%para%group 464 me=grids(auxbas_grid)%pw%pw_grid%para%my_pos 465 glb=rho%pw_grid%bounds(1,:) 466 gub=rho%pw_grid%bounds(2,:) 467 IF ((pos_of_x(glb(1)) .EQ. me).AND.(pos_of_x(gub(1)) .EQ. me)) THEN 468 DO k=lb(3),ub(3) 469 DO j=lb(2),ub(2) 470 grids(auxbas_grid) % pw %cr3d(ub(1)+1,j,k)=rho%cr3d(lb(1),j,k) 471 END DO 472 grids(auxbas_grid) % pw %cr3d(ub(1)+1,ub(2)+1,k)=rho%cr3d(lb(1),lb(2),k) 473 END DO 474 DO j=lb(2),ub(2) 475 grids(auxbas_grid) % pw %cr3d(ub(1)+1,j,ub(3)+1)=rho%cr3d(lb(1),j,lb(3)) 476 END DO 477 grids(auxbas_grid) % pw %cr3d(ub(1)+1,ub(2)+1,ub(3)+1)=rho%cr3d(lb(1),lb(2),lb(3)) 478 ELSE IF (pos_of_x(glb(1)) .EQ. me) THEN 479 ALLOCATE(tmp(rho%pw_grid%bounds_local(1,2):rho%pw_grid%bounds_local(2,2),& 480 rho%pw_grid%bounds_local(1,3):rho%pw_grid%bounds_local(2,3)),& 481 stat=stat) 482 CPASSERT(stat==0) 483 tmp=rho%cr3d(lb(1),:,:) 484 CALL mp_isend(msgin=tmp,dest=pos_of_x(rho%pw_grid%bounds(2,1)),comm=group,& 485 request=request,tag=112) 486 CALL mp_wait(request) 487 ELSE IF (pos_of_x(gub(1)) .EQ. me) THEN 488 ALLOCATE(tmp(rho%pw_grid%bounds_local(1,2):rho%pw_grid%bounds_local(2,2),& 489 rho%pw_grid%bounds_local(1,3):rho%pw_grid%bounds_local(2,3)),& 490 stat=stat) 491 CPASSERT(stat==0) 492 CALL mp_irecv(msgout=tmp,source=pos_of_x(rho%pw_grid%bounds(1,1)),& 493 comm=group,request=request,tag=112) 494 CALL mp_wait(request) 495 496 DO k=lb(3),ub(3) 497 DO j=lb(2),ub(2) 498 grids(auxbas_grid) % pw %cr3d(ub(1)+1,j,k)=tmp(j,k) 499 END DO 500 grids(auxbas_grid) % pw %cr3d(ub(1)+1,ub(2)+1,k)=tmp(lb(2),k) 501 END DO 502 DO j=lb(2),ub(2) 503 grids(auxbas_grid) % pw %cr3d(ub(1)+1,j,ub(3)+1)=tmp(j,lb(3)) 504 END DO 505 grids(auxbas_grid) % pw %cr3d(ub(1)+1,ub(2)+1,ub(3)+1)=tmp(lb(2),lb(3)) 506 END IF 507 IF (ASSOCIATED(tmp)) THEN 508 DEALLOCATE(tmp) 509 END IF 510 ! Further setup of parallelization scheme 511 IF (qmmm_env%par_scheme==do_par_atom) THEN 512 CALL mp_sum(grids(auxbas_grid)%pw%cr3d,para_env%group) 513 END IF 514 ! RealSpace Interpolation 515 CALL section_vals_val_get(interp_section,"kind", i_val=kind_interp) 516 SELECT CASE(kind_interp) 517 CASE(spline3_nopbc_interp, spline3_pbc_interp) 518 ! Spline Interpolator 519 DO Igrid = auxbas_grid, SIZE(grids)-1 520 CALL pw_restrict_s3(grids(Igrid )%pw,& 521 grids(Igrid+1)%pw,& 522 aug_pools(Igrid+1)%pool,& 523 param_section=interp_section) 524 END DO 525 CASE DEFAULT 526 CPABORT("") 527 END SELECT 528 529 shells = .FALSE. 530 CALL qmmm_force_with_gaussian_low(grids, mm_particles,& 531 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index,& 532 qmmm_env%num_mm_atoms, cube_info, para_env, eps_mm_rspace, auxbas_grid, & 533 coarser_grid, qmmm_env%pgfs, qmmm_env%potentials, Forces, aug_pools, & 534 mm_cell, qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%per_potentials, & 535 iw, qmmm_env%par_scheme, qmmm_env%spherical_cutoff, shells) 536 537 IF (qmmm_env%move_mm_charges.OR.qmmm_env%add_mm_charges) THEN 538 CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_charges%added_particles,& 539 qmmm_env%added_charges%mm_atom_chrg,& 540 qmmm_env%added_charges%mm_atom_index, qmmm_env%added_charges%num_mm_atoms,& 541 cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_charges%pgfs,& 542 qmmm_env%added_charges%potentials, Forces_added_charges, aug_pools, mm_cell, & 543 qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%per_potentials, iw, qmmm_env%par_scheme,& 544 qmmm_env%spherical_cutoff, shells) 545 END IF 546 547 IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN 548 shells = .TRUE. 549 CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_shells%added_particles,& 550 qmmm_env%added_shells%mm_core_chrg, & 551 qmmm_env%added_shells%mm_core_index, qmmm_env%added_shells%num_mm_atoms,& 552 cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_shells%pgfs,& 553 qmmm_env%added_shells%potentials, Forces_added_shells, aug_pools, mm_cell, & 554 qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_shells%per_potentials, iw, qmmm_env%par_scheme,& 555 qmmm_env%spherical_cutoff, shells) 556 END IF 557 558 CALL pw_pools_give_back_pws(aug_pools,grids) 559 CALL timestop(handle) 560 561 END SUBROUTINE qmmm_forces_with_gaussian 562 563! ************************************************************************************************** 564!> \brief Evaluates the contribution to the forces due to the 565!> QM/MM potential computed collocating the Electrostatic 566!> Gaussian Potential. Low Level 567!> \param grids ... 568!> \param mm_particles ... 569!> \param mm_charges ... 570!> \param mm_atom_index ... 571!> \param num_mm_atoms ... 572!> \param cube_info ... 573!> \param para_env ... 574!> \param eps_mm_rspace ... 575!> \param auxbas_grid ... 576!> \param coarser_grid ... 577!> \param pgfs ... 578!> \param potentials ... 579!> \param Forces ... 580!> \param aug_pools ... 581!> \param mm_cell ... 582!> \param dOmmOqm ... 583!> \param periodic ... 584!> \param per_potentials ... 585!> \param iw ... 586!> \param par_scheme ... 587!> \param qmmm_spherical_cutoff ... 588!> \param shells ... 589!> \par History 590!> 06.2004 created [tlaino] 591!> \author Teodoro Laino 592! ************************************************************************************************** 593 SUBROUTINE qmmm_force_with_gaussian_low(grids, mm_particles, mm_charges, & 594 mm_atom_index, num_mm_atoms, cube_info, para_env, & 595 eps_mm_rspace, auxbas_grid, coarser_grid, pgfs, potentials, Forces, & 596 aug_pools, mm_cell, dOmmOqm, periodic, per_potentials, iw, par_scheme,& 597 qmmm_spherical_cutoff, shells) 598 TYPE(pw_p_type), DIMENSION(:), POINTER :: grids 599 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 600 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 601 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 602 INTEGER, INTENT(IN) :: num_mm_atoms 603 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info 604 TYPE(cp_para_env_type), POINTER :: para_env 605 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace 606 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid 607 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 608 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: Potentials 609 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces 610 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 611 TYPE(cell_type), POINTER :: mm_cell 612 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 613 LOGICAL, INTENT(in) :: periodic 614 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials 615 INTEGER, INTENT(IN) :: iw, par_scheme 616 REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2) 617 LOGICAL, INTENT(in) :: shells 618 619 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_force_with_gaussian_low', & 620 routineNb = 'qmmm_forces_gaussian_low', routineP = moduleN//':'//routineN 621 622 INTEGER :: handle, handle2, IGauss, ilevel, Imm, & 623 IndMM, IRadTyp, LIndMM, myind, & 624 n_rep_real(3) 625 INTEGER, DIMENSION(2, 3) :: bo 626 REAL(KIND=dp) :: alpha, dvol, height, sph_chrg_factor, W 627 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xdat, ydat, zdat 628 REAL(KIND=dp), DIMENSION(3) :: force, ra 629 TYPE(qmmm_gaussian_type), POINTER :: pgf 630 TYPE(qmmm_per_pot_type), POINTER :: per_pot 631 TYPE(qmmm_pot_type), POINTER :: pot 632 633 CALL timeset(routineN,handle) 634 CALL timeset(routineNb//"_G",handle2) 635 NULLIFY(pgf, pot, per_pot) 636 IF (par_scheme==do_par_atom) myind=0 637 Radius: DO IRadTyp = 1, SIZE(pgfs) 638 pgf => pgfs(IRadTyp)%pgf 639 pot => potentials(IRadTyp)%pot 640 n_rep_real = 0 641 IF (periodic) THEN 642 per_pot => per_potentials(IRadTyp)%pot 643 n_rep_real = per_pot%n_rep_real 644 END IF 645 Gaussian: DO IGauss = 1, pgf%Number_of_Gaussians 646 alpha = 1.0_dp / pgf%Gk(IGauss) 647 alpha = alpha * alpha 648 height = pgf%Ak(IGauss) 649 ilevel = pgf%grid_level(IGauss) 650 dvol = grids(ilevel)%pw%pw_grid%dvol 651 bo = grids(ilevel) % pw % pw_grid % bounds_local 652 ALLOCATE (xdat(2,bo(1,1):bo(2,1))) 653 ALLOCATE (ydat(2,bo(1,2):bo(2,2))) 654 ALLOCATE (zdat(2,bo(1,3):bo(2,3))) 655 Atoms: DO Imm = 1, SIZE(pot%mm_atom_index) 656 IF (par_scheme==do_par_atom) THEN 657 myind = myind + 1 658 IF (MOD(myind,para_env%num_pe)/=para_env%mepos) CYCLE Atoms 659 END IF 660 LIndMM = pot%mm_atom_index(Imm) 661 IndMM = mm_atom_index(LIndMM) 662 IF (shells) THEN 663 ra(:) = pbc(mm_particles(Imm)%r-dOmmOqm, mm_cell)+dOmmOqm 664 ELSE 665 ra(:) = pbc(mm_particles(IndMM)%r-dOmmOqm,mm_cell)+dOmmOqm 666 END IF 667 W = mm_charges(LIndMM) * height 668 force = 0.0_dp 669 ! Possible Spherical Cutoff 670 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 671 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor) 672 W = W * sph_chrg_factor 673 END IF 674 IF (ABS(W)<= EPSILON(0.0_dp)) CYCLE Atoms 675 CALL integrate_gf_rspace_NoPBC(zetp=alpha,& 676 rp=ra,& 677 scale=-1.0_dp,& 678 W=W,& 679 pwgrid=grids(ilevel)%pw,& 680 cube_info=cube_info(ilevel),& 681 eps_mm_rspace=eps_mm_rspace,& 682 xdat=xdat,& 683 ydat=ydat,& 684 zdat=zdat,& 685 bo=bo,& 686 force=force,& 687 n_rep_real=n_rep_real,& 688 mm_cell=mm_cell) 689 force = force * dvol 690 Forces(:,LIndMM) = Forces(:,LIndMM) + force(:) 691 ! 692 ! Debug Statement 693 ! 694 IF (debug_this_module) THEN 695 CALL debug_integrate_gf_rspace_NoPBC(ilevel=ilevel,& 696 zetp=alpha,& 697 rp=ra,& 698 W=W,& 699 pwgrid=grids(ilevel)%pw,& 700 cube_info=cube_info(ilevel),& 701 eps_mm_rspace=eps_mm_rspace,& 702 aug_pools=aug_pools,& 703 debug_force=force,& 704 mm_cell=mm_cell,& 705 auxbas_grid=auxbas_grid,& 706 n_rep_real=n_rep_real,& 707 iw=iw) 708 END IF 709 END DO Atoms 710 DEALLOCATE (xdat) 711 DEALLOCATE (ydat) 712 DEALLOCATE (zdat) 713 END DO Gaussian 714 END DO Radius 715 CALL timestop(handle2) 716 CALL timeset(routineNb//"_R",handle2) 717 IF (periodic) THEN 718 CALL qmmm_forces_with_gaussian_LG (pgfs=pgfs,& 719 cgrid=grids(coarser_grid)%pw,& 720 num_mm_atoms=num_mm_atoms,& 721 mm_charges=mm_charges,& 722 mm_atom_index=mm_atom_index,& 723 mm_particles=mm_particles,& 724 para_env=para_env,& 725 coarser_grid_level=coarser_grid,& 726 Forces=Forces,& 727 per_potentials=per_potentials,& 728 aug_pools=aug_pools,& 729 mm_cell=mm_cell,& 730 dOmmOqm=dOmmOqm,& 731 iw=iw,& 732 par_scheme=par_scheme,& 733 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 734 shells=shells) 735 ELSE 736 CALL qmmm_forces_with_gaussian_LR (pgfs=pgfs,& 737 cgrid=grids(coarser_grid)%pw,& 738 num_mm_atoms=num_mm_atoms,& 739 mm_charges=mm_charges,& 740 mm_atom_index=mm_atom_index,& 741 mm_particles=mm_particles,& 742 para_env=para_env,& 743 coarser_grid_level=coarser_grid,& 744 Forces=Forces,& 745 potentials=potentials,& 746 aug_pools=aug_pools,& 747 mm_cell=mm_cell,& 748 dOmmOqm=dOmmOqm,& 749 iw=iw,& 750 par_scheme=par_scheme,& 751 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 752 shells=shells) 753 END IF 754 CALL timestop(handle2) 755 CALL timestop(handle) 756 END SUBROUTINE qmmm_force_with_gaussian_low 757 758! ************************************************************************************************** 759!> \brief Evaluates the contribution to the forces due to the Long Range 760!> part of the QM/MM potential computed collocating the Electrostatic 761!> Gaussian Potential. 762!> \param pgfs ... 763!> \param cgrid ... 764!> \param num_mm_atoms ... 765!> \param mm_charges ... 766!> \param mm_atom_index ... 767!> \param mm_particles ... 768!> \param para_env ... 769!> \param coarser_grid_level ... 770!> \param Forces ... 771!> \param per_potentials ... 772!> \param aug_pools ... 773!> \param mm_cell ... 774!> \param dOmmOqm ... 775!> \param iw ... 776!> \param par_scheme ... 777!> \param qmmm_spherical_cutoff ... 778!> \param shells ... 779!> \par History 780!> 08.2004 created [tlaino] 781!> \author Teodoro Laino 782! ************************************************************************************************** 783 SUBROUTINE qmmm_forces_with_gaussian_LG (pgfs,cgrid,num_mm_atoms,mm_charges,mm_atom_index,& 784 mm_particles,para_env, coarser_grid_level,Forces, per_potentials,& 785 aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells) 786 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 787 TYPE(pw_type), POINTER :: cgrid 788 INTEGER, INTENT(IN) :: num_mm_atoms 789 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 790 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 791 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 792 TYPE(cp_para_env_type), POINTER :: para_env 793 INTEGER, INTENT(IN) :: coarser_grid_level 794 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces 795 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials 796 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 797 TYPE(cell_type), POINTER :: mm_cell 798 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 799 INTEGER, INTENT(IN) :: iw, par_scheme 800 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff 801 LOGICAL :: shells 802 803 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces_with_gaussian_LG', & 804 routineP = moduleN//':'//routineN 805 806 INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, ij3, ij4, ik1, ik2, ik3, ik4, Imm, & 807 IndMM, IRadTyp, ivec(3), j, k, LIndMM, my_i, my_j, my_k, myind, npts(3) 808 INTEGER, DIMENSION(2, 3) :: bo, gbo 809 REAL(KIND=dp) :: a1, a2, a3, abc_X(4,4), abc_X_Y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, & 810 dr1, dr1c, dr1i, dr2, dr2c, dr2i, dr3, dr3c, dr3i, dvol, e1, e2, e3, f1, f2, f3, fac, & 811 ft1, ft2, ft3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, & 812 rt3, rv1, rv2, rv3, s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, & 813 sph_chrg_factor, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, & 814 v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, v4o, xd1, xd2, xd3, xs1, xs2, xs3 815 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: LForces 816 REAL(KIND=dp), DIMENSION(3) :: ra, val, vec 817 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid, grid2 818 TYPE(pw_type), POINTER :: pw 819 TYPE(qmmm_per_pot_type), POINTER :: per_pot 820 821 CALL timeset(routineN,handle) 822 NULLIFY(grid) 823 ALLOCATE(LForces(3,num_mm_atoms)) 824 LForces = 0.0_dp 825 dr1c = cgrid%pw_grid%dr(1) 826 dr2c = cgrid%pw_grid%dr(2) 827 dr3c = cgrid%pw_grid%dr(3) 828 dvol = cgrid%pw_grid%dvol 829 gbo = cgrid%pw_grid%bounds 830 bo = cgrid%pw_grid%bounds_local 831 grid => cgrid%cr3d 832 IF (par_scheme==do_par_atom) myind = 0 833 Radius: DO IRadTyp = 1, SIZE(pgfs) 834 per_pot => per_potentials(IRadTyp)%pot 835 pw => per_pot%TabLR 836 grid2 => pw%cr3d(:,:,:) 837 npts = pw%pw_grid%npts 838 dr1 = pw%pw_grid%dr(1) 839 dr2 = pw%pw_grid%dr(2) 840 dr3 = pw%pw_grid%dr(3) 841 dr1i = 1.0_dp / dr1 842 dr2i = 1.0_dp / dr2 843 dr3i = 1.0_dp / dr3 844 Atoms: DO Imm = 1, SIZE(per_pot%mm_atom_index) 845 IF (par_scheme==do_par_atom) THEN 846 myind = myind + 1 847 IF (MOD(myind,para_env%num_pe)/=para_env%mepos) CYCLE 848 END IF 849 LIndMM = per_pot%mm_atom_index(Imm) 850 IndMM = mm_atom_index(LIndMM) 851 IF (shells) THEN 852 ra(:) = pbc(mm_particles(LIndMM)%r-dOmmOqm, mm_cell)+dOmmOqm 853 ELSE 854 ra(:) = pbc(mm_particles(IndMM)%r-dOmmOqm,mm_cell)+dOmmOqm 855 END IF 856 qt = mm_charges(LIndMM) 857 ! Possible Spherical Cutoff 858 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 859 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor) 860 qt = qt * sph_chrg_factor 861 END IF 862 IF (ABS(qt)<= EPSILON(0.0_dp)) CYCLE Atoms 863 rt1 = ra(1) 864 rt2 = ra(2) 865 rt3 = ra(3) 866 ft1 = 0.0_dp 867 ft2 = 0.0_dp 868 ft3 = 0.0_dp 869 LoopOnGrid: DO k = bo(1,3), bo(2,3) 870 my_k = k-gbo(1,3) 871 xs3 = REAL(my_k,dp)*dr3c 872 my_j = bo(1,2)-gbo(1,2) 873 xs2 = REAL(my_j,dp)*dr2c 874 rv3 = rt3 - xs3 875 vec(3) = rv3 876 ivec(3) = FLOOR(vec(3)/pw%pw_grid%dr(3)) 877 ik1 = MODULO(ivec(3)-1,npts(3))+1 878 ik2 = MODULO(ivec(3) ,npts(3))+1 879 ik3 = MODULO(ivec(3)+1,npts(3))+1 880 ik4 = MODULO(ivec(3)+2,npts(3))+1 881 xd3 = (vec(3)/dr3)-REAL(ivec(3),kind=dp) 882 p1 = 3.0_dp + xd3 883 p2 = p1*p1 884 p3 = p2*p1 885 q1 = 2.0_dp + xd3 886 q2 = q1*q1 887 q3 = q2*q1 888 r1 = 1.0_dp + xd3 889 r2 = r1*r1 890 r3 = r2*r1 891 u1 = xd3 892 u2 = u1*u1 893 u3 = u2*u1 894 v1o = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3) 895 v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3 896 v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3 897 v4o = 1.0_dp/6.0_dp*u3 898 v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2 899 v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2 900 v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2 901 v4d = 0.5_dp*u2 902 DO j = bo(1,2), bo(2,2) 903 my_i= bo(1,1)-gbo(1,1) 904 xs1 = REAL(my_i,dp)*dr1c 905 rv2 = rt2 - xs2 906 vec(2) = rv2 907 ivec(2) = FLOOR(vec(2)/pw%pw_grid%dr(2)) 908 ij1 = MODULO(ivec(2)-1,npts(2))+1 909 ij2 = MODULO(ivec(2) ,npts(2))+1 910 ij3 = MODULO(ivec(2)+1,npts(2))+1 911 ij4 = MODULO(ivec(2)+2,npts(2))+1 912 xd2 = (vec(2)/dr2)-REAL(ivec(2),kind=dp) 913 e1 = 3.0_dp + xd2 914 e2 = e1*e1 915 e3 = e2*e1 916 f1 = 2.0_dp + xd2 917 f2 = f1*f1 918 f3 = f2*f1 919 g1 = 1.0_dp + xd2 920 g2 = g1*g1 921 g3 = g2*g1 922 h1 = xd2 923 h2 = h1*h1 924 h3 = h2*h1 925 s1o = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3) 926 s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3 927 s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3 928 s4o = 1.0_dp/6.0_dp*h3 929 s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2 930 s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2 931 s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2 932 s4d = 0.5_dp*h2 933 DO i = bo(1,1), bo(2,1) 934 rv1 = rt1 - xs1 935 vec(1) = rv1 936 ivec(1) = FLOOR(vec(1)/pw%pw_grid%dr(1)) 937 ii1 = MODULO(ivec(1)-1,npts(1))+1 938 ii2 = MODULO(ivec(1) ,npts(1))+1 939 ii3 = MODULO(ivec(1)+1,npts(1))+1 940 ii4 = MODULO(ivec(1)+2,npts(1))+1 941 xd1 = (vec(1)/dr1)-REAL(ivec(1),kind=dp) 942 a1 = 3.0_dp + xd1 943 a2 = a1*a1 944 a3 = a2*a1 945 b1 = 2.0_dp + xd1 946 b2 = b1*b1 947 b3 = b2*b1 948 c1 = 1.0_dp + xd1 949 c2 = c1*c1 950 c3 = c2*c1 951 d1 = xd1 952 d2 = d1*d1 953 d3 = d2*d1 954 t1o = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3) 955 t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3 956 t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3 957 t4o = 1.0_dp/6.0_dp*d3 958 t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2 959 t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2 960 t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2 961 t4d = 0.5_dp*d2 962 963 964!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 965! 966! # v then t then s 967! 968! for ii in 1 2 3; do 969! if [[ $ii == 1 ]]; then ld=t; fi 970! if [[ $ii == 2 ]]; then ld=s; fi 971! if [[ $ii == 3 ]]; then ld=v; fi 972! # 973! for l in t s v; do 974! for i in 1 2 3 4; do 975! if [[ $ld == $l ]]; then 976! echo "$l$i = $l${i}d*dr${ii}i" 977! else 978! echo "$l$i = $l${i}o" 979! fi 980! done 981! done 982! # 983! for i in 1 2 3 4; do 984! for j in 1 2 3 4; do 985! echo -n "abc_X($i,$j) = " 986! for k in 1 2 3 4; do 987! if [ $k == 4 ]; then 988! echo "grid2(ii$i,ij$j,ik$k)*v$k" 989! else 990! echo -n "grid2(ii$i,ij$j,ik$k)*v$k + " 991! fi 992! done 993! done 994! done 995! echo "" 996! for j in 1 2 3 4; do 997! echo -n "abc_X_Y($j) = " 998! for k in 1 2 3 4; do 999! if [ $k == 4 ]; then 1000! echo "abc_X($k,$j)*t$k" 1001! else 1002! echo -n "abc_X($k,$j)*t$k + " 1003! fi 1004! done 1005! done 1006! echo "" 1007! echo -n "val($ii) = " 1008! for k in 1 2 3 4; do 1009! if [ $k == 4 ]; then 1010! echo "abc_X_Y($k)*s$k" 1011! else 1012! echo -n "abc_X_Y($k)*s$k + " 1013! fi 1014! done 1015! echo "" 1016! done 1017 1018 t1 = t1d*dr1i 1019 t2 = t2d*dr1i 1020 t3 = t3d*dr1i 1021 t4 = t4d*dr1i 1022 s1 = s1o 1023 s2 = s2o 1024 s3 = s3o 1025 s4 = s4o 1026 v1 = v1o 1027 v2 = v2o 1028 v3 = v3o 1029 v4 = v4o 1030 abc_X(1,1) = grid2(ii1,ij1,ik1)*v1 + grid2(ii1,ij1,ik2)*v2 + grid2(ii1,ij1,ik3)*v3 + grid2(ii1,ij1,ik4)*v4 1031 abc_X(1,2) = grid2(ii1,ij2,ik1)*v1 + grid2(ii1,ij2,ik2)*v2 + grid2(ii1,ij2,ik3)*v3 + grid2(ii1,ij2,ik4)*v4 1032 abc_X(1,3) = grid2(ii1,ij3,ik1)*v1 + grid2(ii1,ij3,ik2)*v2 + grid2(ii1,ij3,ik3)*v3 + grid2(ii1,ij3,ik4)*v4 1033 abc_X(1,4) = grid2(ii1,ij4,ik1)*v1 + grid2(ii1,ij4,ik2)*v2 + grid2(ii1,ij4,ik3)*v3 + grid2(ii1,ij4,ik4)*v4 1034 abc_X(2,1) = grid2(ii2,ij1,ik1)*v1 + grid2(ii2,ij1,ik2)*v2 + grid2(ii2,ij1,ik3)*v3 + grid2(ii2,ij1,ik4)*v4 1035 abc_X(2,2) = grid2(ii2,ij2,ik1)*v1 + grid2(ii2,ij2,ik2)*v2 + grid2(ii2,ij2,ik3)*v3 + grid2(ii2,ij2,ik4)*v4 1036 abc_X(2,3) = grid2(ii2,ij3,ik1)*v1 + grid2(ii2,ij3,ik2)*v2 + grid2(ii2,ij3,ik3)*v3 + grid2(ii2,ij3,ik4)*v4 1037 abc_X(2,4) = grid2(ii2,ij4,ik1)*v1 + grid2(ii2,ij4,ik2)*v2 + grid2(ii2,ij4,ik3)*v3 + grid2(ii2,ij4,ik4)*v4 1038 abc_X(3,1) = grid2(ii3,ij1,ik1)*v1 + grid2(ii3,ij1,ik2)*v2 + grid2(ii3,ij1,ik3)*v3 + grid2(ii3,ij1,ik4)*v4 1039 abc_X(3,2) = grid2(ii3,ij2,ik1)*v1 + grid2(ii3,ij2,ik2)*v2 + grid2(ii3,ij2,ik3)*v3 + grid2(ii3,ij2,ik4)*v4 1040 abc_X(3,3) = grid2(ii3,ij3,ik1)*v1 + grid2(ii3,ij3,ik2)*v2 + grid2(ii3,ij3,ik3)*v3 + grid2(ii3,ij3,ik4)*v4 1041 abc_X(3,4) = grid2(ii3,ij4,ik1)*v1 + grid2(ii3,ij4,ik2)*v2 + grid2(ii3,ij4,ik3)*v3 + grid2(ii3,ij4,ik4)*v4 1042 abc_X(4,1) = grid2(ii4,ij1,ik1)*v1 + grid2(ii4,ij1,ik2)*v2 + grid2(ii4,ij1,ik3)*v3 + grid2(ii4,ij1,ik4)*v4 1043 abc_X(4,2) = grid2(ii4,ij2,ik1)*v1 + grid2(ii4,ij2,ik2)*v2 + grid2(ii4,ij2,ik3)*v3 + grid2(ii4,ij2,ik4)*v4 1044 abc_X(4,3) = grid2(ii4,ij3,ik1)*v1 + grid2(ii4,ij3,ik2)*v2 + grid2(ii4,ij3,ik3)*v3 + grid2(ii4,ij3,ik4)*v4 1045 abc_X(4,4) = grid2(ii4,ij4,ik1)*v1 + grid2(ii4,ij4,ik2)*v2 + grid2(ii4,ij4,ik3)*v3 + grid2(ii4,ij4,ik4)*v4 1046 1047 abc_X_Y(1) = abc_X(1,1)*t1 + abc_X(2,1)*t2 + abc_X(3,1)*t3 + abc_X(4,1)*t4 1048 abc_X_Y(2) = abc_X(1,2)*t1 + abc_X(2,2)*t2 + abc_X(3,2)*t3 + abc_X(4,2)*t4 1049 abc_X_Y(3) = abc_X(1,3)*t1 + abc_X(2,3)*t2 + abc_X(3,3)*t3 + abc_X(4,3)*t4 1050 abc_X_Y(4) = abc_X(1,4)*t1 + abc_X(2,4)*t2 + abc_X(3,4)*t3 + abc_X(4,4)*t4 1051 1052 val(1) = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4 1053 1054 t1 = t1o 1055 t2 = t2o 1056 t3 = t3o 1057 t4 = t4o 1058 s1 = s1d*dr2i 1059 s2 = s2d*dr2i 1060 s3 = s3d*dr2i 1061 s4 = s4d*dr2i 1062 !! v1 = v1o 1063 !! v2 = v2o 1064 !! v3 = v3o 1065 !! v4 = v4o 1066 !! abc_X(1,1) = grid2(ii1,ij1,ik1)*v1 + grid2(ii1,ij1,ik2)*v2 + grid2(ii1,ij1,ik3)*v3 + grid2(ii1,ij1,ik4)*v4 1067 !! abc_X(1,2) = grid2(ii1,ij2,ik1)*v1 + grid2(ii1,ij2,ik2)*v2 + grid2(ii1,ij2,ik3)*v3 + grid2(ii1,ij2,ik4)*v4 1068 !! abc_X(1,3) = grid2(ii1,ij3,ik1)*v1 + grid2(ii1,ij3,ik2)*v2 + grid2(ii1,ij3,ik3)*v3 + grid2(ii1,ij3,ik4)*v4 1069 !! abc_X(1,4) = grid2(ii1,ij4,ik1)*v1 + grid2(ii1,ij4,ik2)*v2 + grid2(ii1,ij4,ik3)*v3 + grid2(ii1,ij4,ik4)*v4 1070 !! abc_X(2,1) = grid2(ii2,ij1,ik1)*v1 + grid2(ii2,ij1,ik2)*v2 + grid2(ii2,ij1,ik3)*v3 + grid2(ii2,ij1,ik4)*v4 1071 !! abc_X(2,2) = grid2(ii2,ij2,ik1)*v1 + grid2(ii2,ij2,ik2)*v2 + grid2(ii2,ij2,ik3)*v3 + grid2(ii2,ij2,ik4)*v4 1072 !! abc_X(2,3) = grid2(ii2,ij3,ik1)*v1 + grid2(ii2,ij3,ik2)*v2 + grid2(ii2,ij3,ik3)*v3 + grid2(ii2,ij3,ik4)*v4 1073 !! abc_X(2,4) = grid2(ii2,ij4,ik1)*v1 + grid2(ii2,ij4,ik2)*v2 + grid2(ii2,ij4,ik3)*v3 + grid2(ii2,ij4,ik4)*v4 1074 !! abc_X(3,1) = grid2(ii3,ij1,ik1)*v1 + grid2(ii3,ij1,ik2)*v2 + grid2(ii3,ij1,ik3)*v3 + grid2(ii3,ij1,ik4)*v4 1075 !! abc_X(3,2) = grid2(ii3,ij2,ik1)*v1 + grid2(ii3,ij2,ik2)*v2 + grid2(ii3,ij2,ik3)*v3 + grid2(ii3,ij2,ik4)*v4 1076 !! abc_X(3,3) = grid2(ii3,ij3,ik1)*v1 + grid2(ii3,ij3,ik2)*v2 + grid2(ii3,ij3,ik3)*v3 + grid2(ii3,ij3,ik4)*v4 1077 !! abc_X(3,4) = grid2(ii3,ij4,ik1)*v1 + grid2(ii3,ij4,ik2)*v2 + grid2(ii3,ij4,ik3)*v3 + grid2(ii3,ij4,ik4)*v4 1078 !! abc_X(4,1) = grid2(ii4,ij1,ik1)*v1 + grid2(ii4,ij1,ik2)*v2 + grid2(ii4,ij1,ik3)*v3 + grid2(ii4,ij1,ik4)*v4 1079 !! abc_X(4,2) = grid2(ii4,ij2,ik1)*v1 + grid2(ii4,ij2,ik2)*v2 + grid2(ii4,ij2,ik3)*v3 + grid2(ii4,ij2,ik4)*v4 1080 !! abc_X(4,3) = grid2(ii4,ij3,ik1)*v1 + grid2(ii4,ij3,ik2)*v2 + grid2(ii4,ij3,ik3)*v3 + grid2(ii4,ij3,ik4)*v4 1081 !! abc_X(4,4) = grid2(ii4,ij4,ik1)*v1 + grid2(ii4,ij4,ik2)*v2 + grid2(ii4,ij4,ik3)*v3 + grid2(ii4,ij4,ik4)*v4 1082 1083 abc_X_Y(1) = abc_X(1,1)*t1 + abc_X(2,1)*t2 + abc_X(3,1)*t3 + abc_X(4,1)*t4 1084 abc_X_Y(2) = abc_X(1,2)*t1 + abc_X(2,2)*t2 + abc_X(3,2)*t3 + abc_X(4,2)*t4 1085 abc_X_Y(3) = abc_X(1,3)*t1 + abc_X(2,3)*t2 + abc_X(3,3)*t3 + abc_X(4,3)*t4 1086 abc_X_Y(4) = abc_X(1,4)*t1 + abc_X(2,4)*t2 + abc_X(3,4)*t3 + abc_X(4,4)*t4 1087 1088 val(2) = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4 1089 1090 t1 = t1o 1091 t2 = t2o 1092 t3 = t3o 1093 t4 = t4o 1094 s1 = s1o 1095 s2 = s2o 1096 s3 = s3o 1097 s4 = s4o 1098 v1 = v1d*dr3i 1099 v2 = v2d*dr3i 1100 v3 = v3d*dr3i 1101 v4 = v4d*dr3i 1102 abc_X(1,1) = grid2(ii1,ij1,ik1)*v1 + grid2(ii1,ij1,ik2)*v2 + grid2(ii1,ij1,ik3)*v3 + grid2(ii1,ij1,ik4)*v4 1103 abc_X(1,2) = grid2(ii1,ij2,ik1)*v1 + grid2(ii1,ij2,ik2)*v2 + grid2(ii1,ij2,ik3)*v3 + grid2(ii1,ij2,ik4)*v4 1104 abc_X(1,3) = grid2(ii1,ij3,ik1)*v1 + grid2(ii1,ij3,ik2)*v2 + grid2(ii1,ij3,ik3)*v3 + grid2(ii1,ij3,ik4)*v4 1105 abc_X(1,4) = grid2(ii1,ij4,ik1)*v1 + grid2(ii1,ij4,ik2)*v2 + grid2(ii1,ij4,ik3)*v3 + grid2(ii1,ij4,ik4)*v4 1106 abc_X(2,1) = grid2(ii2,ij1,ik1)*v1 + grid2(ii2,ij1,ik2)*v2 + grid2(ii2,ij1,ik3)*v3 + grid2(ii2,ij1,ik4)*v4 1107 abc_X(2,2) = grid2(ii2,ij2,ik1)*v1 + grid2(ii2,ij2,ik2)*v2 + grid2(ii2,ij2,ik3)*v3 + grid2(ii2,ij2,ik4)*v4 1108 abc_X(2,3) = grid2(ii2,ij3,ik1)*v1 + grid2(ii2,ij3,ik2)*v2 + grid2(ii2,ij3,ik3)*v3 + grid2(ii2,ij3,ik4)*v4 1109 abc_X(2,4) = grid2(ii2,ij4,ik1)*v1 + grid2(ii2,ij4,ik2)*v2 + grid2(ii2,ij4,ik3)*v3 + grid2(ii2,ij4,ik4)*v4 1110 abc_X(3,1) = grid2(ii3,ij1,ik1)*v1 + grid2(ii3,ij1,ik2)*v2 + grid2(ii3,ij1,ik3)*v3 + grid2(ii3,ij1,ik4)*v4 1111 abc_X(3,2) = grid2(ii3,ij2,ik1)*v1 + grid2(ii3,ij2,ik2)*v2 + grid2(ii3,ij2,ik3)*v3 + grid2(ii3,ij2,ik4)*v4 1112 abc_X(3,3) = grid2(ii3,ij3,ik1)*v1 + grid2(ii3,ij3,ik2)*v2 + grid2(ii3,ij3,ik3)*v3 + grid2(ii3,ij3,ik4)*v4 1113 abc_X(3,4) = grid2(ii3,ij4,ik1)*v1 + grid2(ii3,ij4,ik2)*v2 + grid2(ii3,ij4,ik3)*v3 + grid2(ii3,ij4,ik4)*v4 1114 abc_X(4,1) = grid2(ii4,ij1,ik1)*v1 + grid2(ii4,ij1,ik2)*v2 + grid2(ii4,ij1,ik3)*v3 + grid2(ii4,ij1,ik4)*v4 1115 abc_X(4,2) = grid2(ii4,ij2,ik1)*v1 + grid2(ii4,ij2,ik2)*v2 + grid2(ii4,ij2,ik3)*v3 + grid2(ii4,ij2,ik4)*v4 1116 abc_X(4,3) = grid2(ii4,ij3,ik1)*v1 + grid2(ii4,ij3,ik2)*v2 + grid2(ii4,ij3,ik3)*v3 + grid2(ii4,ij3,ik4)*v4 1117 abc_X(4,4) = grid2(ii4,ij4,ik1)*v1 + grid2(ii4,ij4,ik2)*v2 + grid2(ii4,ij4,ik3)*v3 + grid2(ii4,ij4,ik4)*v4 1118 1119 abc_X_Y(1) = abc_X(1,1)*t1 + abc_X(2,1)*t2 + abc_X(3,1)*t3 + abc_X(4,1)*t4 1120 abc_X_Y(2) = abc_X(1,2)*t1 + abc_X(2,2)*t2 + abc_X(3,2)*t3 + abc_X(4,2)*t4 1121 abc_X_Y(3) = abc_X(1,3)*t1 + abc_X(2,3)*t2 + abc_X(3,3)*t3 + abc_X(4,3)*t4 1122 abc_X_Y(4) = abc_X(1,4)*t1 + abc_X(2,4)*t2 + abc_X(3,4)*t3 + abc_X(4,4)*t4 1123 1124 val(3) = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4 1125!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1126 1127 fac = grid(i,j,k) 1128 ft1 = ft1 + val(1) * fac 1129 ft2 = ft2 + val(2) * fac 1130 ft3 = ft3 + val(3) * fac 1131 xs1 = xs1 + dr1c 1132 END DO 1133 xs2 = xs2 + dr2c 1134 END DO 1135 END DO LoopOnGrid 1136 qt = - qt * dvol 1137 LForces(1,LindMM) = ft1 * qt 1138 LForces(2,LindMM) = ft2 * qt 1139 LForces(3,LindMM) = ft3 * qt 1140 1141 Forces(1,LIndMM) = Forces(1,LIndMM) + LForces(1,LindMM) 1142 Forces(2,LIndMM) = Forces(2,LIndMM) + LForces(2,LindMM) 1143 Forces(3,LIndMM) = Forces(3,LIndMM) + LForces(3,LindMM) 1144 END DO Atoms 1145 END DO Radius 1146 ! 1147 ! Debug Statement 1148 ! 1149 IF (debug_this_module) THEN 1150 CALL debug_qmmm_forces_with_gauss_LG(pgfs=pgfs,& 1151 aug_pools=aug_pools,& 1152 rho=cgrid,& 1153 num_mm_atoms=num_mm_atoms,& 1154 mm_charges=mm_charges,& 1155 mm_atom_index=mm_atom_index,& 1156 mm_particles=mm_particles,& 1157 coarser_grid_level=coarser_grid_level,& 1158 debug_force=LForces,& 1159 per_potentials=per_potentials,& 1160 para_env=para_env,& 1161 mm_cell=mm_cell,& 1162 dOmmOqm=dOmmOqm,& 1163 iw=iw,& 1164 par_scheme=par_scheme,& 1165 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 1166 shells=shells) 1167 END IF 1168 DEALLOCATE(LForces) 1169 CALL timestop(handle) 1170 END SUBROUTINE qmmm_forces_with_gaussian_LG 1171 1172! ************************************************************************************************** 1173!> \brief Evaluates the contribution to the forces due to the Long Range 1174!> part of the QM/MM potential computed collocating the Electrostatic 1175!> Gaussian Potential. 1176!> \param pgfs ... 1177!> \param cgrid ... 1178!> \param num_mm_atoms ... 1179!> \param mm_charges ... 1180!> \param mm_atom_index ... 1181!> \param mm_particles ... 1182!> \param para_env ... 1183!> \param coarser_grid_level ... 1184!> \param Forces ... 1185!> \param potentials ... 1186!> \param aug_pools ... 1187!> \param mm_cell ... 1188!> \param dOmmOqm ... 1189!> \param iw ... 1190!> \param par_scheme ... 1191!> \param qmmm_spherical_cutoff ... 1192!> \param shells ... 1193!> \par History 1194!> 08.2004 created [tlaino] 1195!> \author Teodoro Laino 1196! ************************************************************************************************** 1197 SUBROUTINE qmmm_forces_with_gaussian_LR (pgfs,cgrid,num_mm_atoms,mm_charges,mm_atom_index,& 1198 mm_particles,para_env, coarser_grid_level,Forces, potentials,& 1199 aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells) 1200 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 1201 TYPE(pw_type), POINTER :: cgrid 1202 INTEGER, INTENT(IN) :: num_mm_atoms 1203 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 1204 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1205 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 1206 TYPE(cp_para_env_type), POINTER :: para_env 1207 INTEGER, INTENT(IN) :: coarser_grid_level 1208 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces 1209 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: Potentials 1210 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 1211 TYPE(cell_type), POINTER :: mm_cell 1212 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 1213 INTEGER, INTENT(IN) :: iw, par_scheme 1214 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff 1215 LOGICAL :: shells 1216 1217 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_forces_with_gaussian_LR', & 1218 routineP = moduleN//':'//routineN 1219 1220 INTEGER :: handle, i, Imm, IndMM, IRadTyp, ix, j, & 1221 k, LIndMM, my_i, my_j, my_k, myind, & 1222 n1, n2, n3 1223 INTEGER, DIMENSION(2, 3) :: bo, gbo 1224 REAL(KIND=dp) :: dr1, dr2, dr3, dvol, dx, fac, ft1, ft2, & 1225 ft3, qt, r, r2, rd1, rd2, rd3, rt1, & 1226 rt2, rt3, rv1, rv2, rv3, rx, rx2, & 1227 sph_chrg_factor, Term, xs1, xs2, xs3 1228 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: LForces 1229 REAL(KIND=dp), DIMENSION(3) :: ra 1230 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pot0_2 1231 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid 1232 TYPE(qmmm_pot_type), POINTER :: pot 1233 1234 CALL timeset(routineN,handle) 1235 ALLOCATE(LForces(3,num_mm_atoms)) 1236 LForces = 0.0_dp 1237 n1 = cgrid%pw_grid%npts(1) 1238 n2 = cgrid%pw_grid%npts(2) 1239 n3 = cgrid%pw_grid%npts(3) 1240 dr1 = cgrid%pw_grid%dr(1) 1241 dr2 = cgrid%pw_grid%dr(2) 1242 dr3 = cgrid%pw_grid%dr(3) 1243 dvol = cgrid%pw_grid%dvol 1244 gbo = cgrid%pw_grid%bounds 1245 bo = cgrid%pw_grid%bounds_local 1246 grid => cgrid%cr3d 1247 IF (par_scheme==do_par_atom) myind = 0 1248 Radius: DO IRadTyp = 1, SIZE(pgfs) 1249 pot => potentials(IRadTyp)%pot 1250 dx = Pot%dx 1251 pot0_2 => Pot%pot0_2 1252 Atoms: DO Imm = 1, SIZE(pot%mm_atom_index) 1253 IF (par_scheme==do_par_atom) THEN 1254 myind = myind + 1 1255 IF (MOD(myind,para_env%num_pe)/=para_env%mepos) CYCLE 1256 END IF 1257 LIndMM = pot%mm_atom_index(Imm) 1258 IndMM = mm_atom_index(LIndMM) 1259 ra(:) = pbc(mm_particles(IndMM)%r-dOmmOqm,mm_cell)+dOmmOqm 1260 IF (shells) & 1261 ra(:) = pbc(mm_particles(LIndMM)%r-dOmmOqm, mm_cell)+dOmmOqm 1262 qt = mm_charges(LIndMM) 1263 ! Possible Spherical Cutoff 1264 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 1265 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor) 1266 qt = qt * sph_chrg_factor 1267 END IF 1268 IF (ABS(qt)<= EPSILON(0.0_dp)) CYCLE Atoms 1269 rt1 = ra(1) 1270 rt2 = ra(2) 1271 rt3 = ra(3) 1272 ft1 = 0.0_dp 1273 ft2 = 0.0_dp 1274 ft3 = 0.0_dp 1275 LoopOnGrid: DO k = bo(1,3), bo(2,3) 1276 my_k = k-gbo(1,3) 1277 xs3 = REAL(my_k,dp)*dr3 1278 my_j = bo(1,2)-gbo(1,2) 1279 xs2 = REAL(my_j,dp)*dr2 1280 rv3 = rt3 - xs3 1281 DO j = bo(1,2), bo(2,2) 1282 my_i= bo(1,1)-gbo(1,1) 1283 xs1 = REAL(my_i,dp)*dr1 1284 rv2 = rt2 - xs2 1285 DO i = bo(1,1), bo(2,1) 1286 rv1 = rt1 - xs1 1287 r2 = rv1*rv1 + rv2*rv2 + rv3*rv3 1288 r = SQRT(r2) 1289 ix = FLOOR(r/dx)+1 1290 rx = (r-REAL(ix-1,dp)*dx)/dx 1291 rx2 = rx*rx 1292 Term = pot0_2(1,ix )*(-6.0_dp*(rx-rx2)) & 1293 +pot0_2(2,ix )*(1.0_dp-4.0_dp*rx+3.0_dp*rx2) & 1294 +pot0_2(1,ix+1)*(6.0_dp*(rx-rx2)) & 1295 +pot0_2(2,ix+1)*(-2.0_dp*rx+3.0_dp*rx2) 1296 fac = grid(i,j,k) * Term 1297 IF ( r == 0.0_dp ) THEN 1298 rd1 = 1.0_dp 1299 rd2 = 1.0_dp 1300 rd3 = 1.0_dp 1301 ELSE 1302 rd1 = rv1 / r 1303 rd2 = rv2 / r 1304 rd3 = rv3 / r 1305 END IF 1306 ft1 = ft1 + fac * rd1 1307 ft2 = ft2 + fac * rd2 1308 ft3 = ft3 + fac * rd3 1309 xs1 = xs1 + dr1 1310 END DO 1311 xs2 = xs2 + dr2 1312 END DO 1313 END DO LoopOnGrid 1314 qt = - qt * dvol / dx 1315 LForces(1,LindMM) = ft1 * qt 1316 LForces(2,LindMM) = ft2 * qt 1317 LForces(3,LindMM) = ft3 * qt 1318 1319 Forces(1,LIndMM) = Forces(1,LIndMM) + LForces(1,LindMM) 1320 Forces(2,LIndMM) = Forces(2,LIndMM) + LForces(2,LindMM) 1321 Forces(3,LIndMM) = Forces(3,LIndMM) + LForces(3,LindMM) 1322 END DO Atoms 1323 END DO Radius 1324 ! 1325 ! Debug Statement 1326 ! 1327 IF (debug_this_module) THEN 1328 CALL debug_qmmm_forces_with_gauss_LR(pgfs=pgfs,& 1329 aug_pools=aug_pools,& 1330 rho=cgrid,& 1331 num_mm_atoms=num_mm_atoms,& 1332 mm_charges=mm_charges,& 1333 mm_atom_index=mm_atom_index,& 1334 mm_particles=mm_particles,& 1335 coarser_grid_level=coarser_grid_level,& 1336 debug_force=LForces,& 1337 potentials=potentials,& 1338 para_env=para_env,& 1339 mm_cell=mm_cell,& 1340 dOmmOqm=dOmmOqm,& 1341 iw=iw,& 1342 par_scheme=par_scheme,& 1343 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 1344 shells=shells) 1345 END IF 1346 1347 DEALLOCATE(LForces) 1348 CALL timestop(handle) 1349 END SUBROUTINE qmmm_forces_with_gaussian_LR 1350 1351! ************************************************************************************************** 1352!> \brief Evaluates numerically QM/MM forces and compares them with 1353!> the analytically computed ones. 1354!> It is evaluated only when debug_this_module is set to .TRUE. 1355!> \param rho ... 1356!> \param qs_env ... 1357!> \param qmmm_env ... 1358!> \param Analytical_Forces ... 1359!> \param mm_particles ... 1360!> \param mm_atom_index ... 1361!> \param num_mm_atoms ... 1362!> \param interp_section ... 1363!> \param mm_cell ... 1364!> \par History 1365!> 08.2004 created [tlaino] 1366!> \author Teodoro Laino 1367! ************************************************************************************************** 1368 SUBROUTINE qmmm_debug_forces(rho,qs_env,qmmm_env,Analytical_Forces,& 1369 mm_particles,mm_atom_index,num_mm_atoms,& 1370 interp_section,mm_cell) 1371 TYPE(pw_type), POINTER :: rho 1372 TYPE(qs_environment_type), POINTER :: qs_env 1373 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 1374 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Analytical_Forces 1375 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 1376 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1377 INTEGER, INTENT(IN) :: num_mm_atoms 1378 TYPE(section_vals_type), POINTER :: interp_section 1379 TYPE(cell_type), POINTER :: mm_cell 1380 1381 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_debug_forces', & 1382 routineP = moduleN//':'//routineN 1383 1384 INTEGER :: handle, I, IndMM, iw, J, K 1385 REAL(KIND=dp) :: Coord_save 1386 REAL(KIND=dp), DIMENSION(2) :: energy 1387 REAL(KIND=dp), DIMENSION(3) :: Err 1388 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Num_Forces 1389 TYPE(cp_logger_type), POINTER :: logger 1390 TYPE(cp_para_env_type), POINTER :: para_env 1391 TYPE(pw_env_type), POINTER :: pw_env 1392 TYPE(pw_p_type) :: v_qmmm_rspace 1393 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 1394 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc 1395 TYPE(section_vals_type), POINTER :: input_section, print_section 1396 1397 CALL timeset(routineN,handle) 1398 NULLIFY( Num_Forces ) 1399 CALL get_qs_env(qs_env=qs_env,& 1400 pw_env=pw_env,& 1401 input=input_section,& 1402 para_env=para_env) 1403 1404 print_section => section_vals_get_subs_vals(input_section,"QMMM%PRINT") 1405 logger => cp_get_default_logger() 1406 iw=cp_print_key_unit_nr(logger,print_section,"PROGRAM_RUN_INFO",extension=".qmmmLog") 1407 CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools) 1408 CALL pw_pool_create_pw(pw_pools(1)%pool, v_qmmm_rspace%pw,& 1409 use_data=REALDATA3D, in_space=REALSPACE) 1410 ALLOCATE(Num_Forces(3,num_mm_atoms)) 1411 ks_qmmm_env_loc => qs_env%ks_qmmm_env 1412 IF (iw>0) WRITE(iw,'(/A)')"DEBUG SECTION:" 1413 Atoms: DO I = 1, num_mm_atoms 1414 IndMM = mm_atom_index(I) 1415 Coords: DO J = 1, 3 1416 Coord_save = mm_particles(IndMM)%r(J) 1417 energy = 0.0_dp 1418 Diff: DO K = 1, 2 1419 mm_particles(IndMM)%r(J) = Coord_save + (-1)**K * Dx 1420 CALL pw_zero(v_qmmm_rspace%pw) 1421 SELECT CASE(qmmm_env%qmmm_coupl_type) 1422 CASE(do_qmmm_coulomb) 1423 CPABORT("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.") 1424 CASE(do_qmmm_pcharge) 1425 CPABORT("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.") 1426 CASE(do_qmmm_gauss,do_qmmm_swave) 1427 CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env,& 1428 v_qmmm=v_qmmm_rspace,& 1429 mm_particles=mm_particles,& 1430 aug_pools=qmmm_env%aug_pools,& 1431 para_env=para_env,& 1432 eps_mm_rspace=qmmm_env%eps_mm_rspace,& 1433 cube_info=ks_qmmm_env_loc%cube_info,& 1434 pw_pools=pw_pools,& 1435 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid,& 1436 coarser_grid=qmmm_env%gridlevel_info%coarser_grid,& 1437 interp_section=interp_section,& 1438 mm_cell=mm_cell) 1439 CASE(do_qmmm_none) 1440 CYCLE Diff 1441 CASE DEFAULT 1442 IF (iw>0) WRITE(iw,'(T3,A)')"Unknown Coupling..." 1443 CPABORT("") 1444 END SELECT 1445 energy(K) = pw_integral_ab ( rho, v_qmmm_rspace%pw) 1446 END DO Diff 1447 IF (iw>0) THEN 1448 WRITE(iw,'(A,I6,A,I3,A,2F15.9)')& 1449 "DEBUG :: MM Atom = ",IndMM," Coord = ",J," Energies (+/-) :: ",energy(2), energy(1) 1450 END IF 1451 Num_Forces(J,I) = ( energy(2) - energy(1) ) / (2.0_dp * Dx) 1452 mm_particles(IndMM)%r(J) = Coord_save 1453 END DO Coords 1454 END DO Atoms 1455 1456 SELECT CASE(qmmm_env%qmmm_coupl_type) 1457 CASE(do_qmmm_coulomb) 1458 CPABORT("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.") 1459 CASE(do_qmmm_pcharge) 1460 CPABORT("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.") 1461 CASE(do_qmmm_gauss,do_qmmm_swave) 1462 IF (iw>0) WRITE(iw,'(/A/)')"CHECKING NUMERICAL Vs ANALYTICAL FORCES (Err%):" 1463 DO I = 1, num_mm_atoms 1464 IndMM = mm_atom_index(I) 1465 Err = 0.0_dp 1466 DO K = 1, 3 1467 IF (ABS(Num_Forces(K,I)) >= 5.0E-5_dp) THEN 1468 Err(K) = (Analytical_Forces(K,I)-Num_Forces(K,I))/Num_Forces(K,I)*100.0_dp 1469 END IF 1470 END DO 1471 IF (iw>0)& 1472 WRITE(iw,100)IndMM,Analytical_Forces(1,I),Num_Forces(1,I),Err(1),& 1473 Analytical_Forces(2,I),Num_Forces(2,I),Err(2),& 1474 Analytical_Forces(3,I),Num_Forces(3,I),Err(3) 1475 CPASSERT(ABS(Err(1))<=MaxErr) 1476 CPASSERT(ABS(Err(2))<=MaxErr) 1477 CPASSERT(ABS(Err(3))<=MaxErr) 1478 END DO 1479 CASE(do_qmmm_none) 1480 IF (iw>0) WRITE(iw,'(T3,A)')"No QM/MM Derivatives to debug. Just Mechanical Coupling!" 1481 CASE DEFAULT 1482 IF (iw>0) WRITE(iw,'(T3,A)')"Unknown Coupling..." 1483 CPABORT("") 1484 END SELECT 1485 CALL cp_print_key_finished_output(iw,logger,print_section,"PROGRAM_RUN_INFO") 1486 1487 CALL pw_pool_give_back_pw ( pw_pools(1)%pool, v_qmmm_rspace%pw) 1488 DEALLOCATE(Num_Forces) 1489 CALL timestop(handle) 1490100 FORMAT(I5,2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ") 1491 END SUBROUTINE qmmm_debug_forces 1492 1493! ************************************************************************************************** 1494!> \brief Debugs the integrate_gf_rspace_NoPBC.. It may helps ;-P 1495!> \param ilevel ... 1496!> \param zetp ... 1497!> \param rp ... 1498!> \param W ... 1499!> \param pwgrid ... 1500!> \param cube_info ... 1501!> \param eps_mm_rspace ... 1502!> \param aug_pools ... 1503!> \param debug_force ... 1504!> \param mm_cell ... 1505!> \param auxbas_grid ... 1506!> \param n_rep_real ... 1507!> \param iw ... 1508!> \par History 1509!> 08.2004 created [tlaino] 1510!> \author Teodoro Laino 1511! ************************************************************************************************** 1512 SUBROUTINE debug_integrate_gf_rspace_NoPBC(ilevel, zetp, rp, W, pwgrid, cube_info,& 1513 eps_mm_rspace, aug_pools, debug_force,& 1514 mm_cell,auxbas_grid, n_rep_real, iw) 1515 INTEGER, INTENT(IN) :: ilevel 1516 REAL(KIND=dp), INTENT(IN) :: zetp 1517 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rp 1518 REAL(KIND=dp), INTENT(IN) :: W 1519 TYPE(pw_type), POINTER :: pwgrid 1520 TYPE(cube_info_type), INTENT(IN) :: cube_info 1521 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace 1522 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 1523 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: debug_force 1524 TYPE(cell_type), POINTER :: mm_cell 1525 INTEGER, INTENT(IN) :: auxbas_grid 1526 INTEGER, DIMENSION(3), INTENT(IN) :: n_rep_real 1527 INTEGER, INTENT(IN) :: iw 1528 1529 CHARACTER(len=*), PARAMETER :: routineN = 'debug_integrate_gf_rspace_NoPBC', & 1530 routineP = moduleN//':'//routineN 1531 1532 INTEGER :: handle, i, igrid, k, ngrids 1533 INTEGER, DIMENSION(2, 3) :: bo2 1534 INTEGER, SAVE :: Icount 1535 REAL(KIND=dp), DIMENSION(2) :: energy 1536 REAL(KIND=dp), DIMENSION(3) :: Err, force, myrp 1537 REAL(KIND=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat 1538 TYPE(pw_p_type), DIMENSION(:), POINTER :: grids 1539 1540 DATA Icount /0/ 1541 ! Statements 1542 CALL timeset(routineN,handle) 1543 NULLIFY(grids) 1544 !Statements 1545 ngrids = SIZE(aug_pools) 1546 CALL pw_pools_create_pws(aug_pools,grids,use_data=REALDATA3D,in_space=REALSPACE) 1547 DO igrid=1,ngrids 1548 CALL pw_zero(grids(igrid)%pw) 1549 END DO 1550 bo2 = grids(auxbas_grid)%pw%pw_grid%bounds 1551 ALLOCATE (xdat(bo2(1,1):bo2(2,1))) 1552 ALLOCATE (ydat(bo2(1,2):bo2(2,2))) 1553 ALLOCATE (zdat(bo2(1,3):bo2(2,3))) 1554 1555 Icount = Icount + 1 1556 DO i = 1, 3 1557 DO k = 1, 2 1558 myrp = rp 1559 myrp(i) = myrp(i) + (-1.0_dp)**k * Dx 1560 CALL pw_zero(grids(ilevel)%pw) 1561 CALL collocate_gf_rspace_NoPBC(zetp=zetp,& 1562 rp=myrp,& 1563 scale=-1.0_dp,& 1564 W=W,& 1565 pwgrid=grids(ilevel)%pw,& 1566 cube_info=cube_info,& 1567 eps_mm_rspace=eps_mm_rspace,& 1568 xdat=xdat,& 1569 ydat=ydat,& 1570 zdat=zdat,& 1571 bo2=bo2,& 1572 n_rep_real=n_rep_real,& 1573 mm_cell=mm_cell) 1574 1575 energy(k) = pw_integral_ab(pwgrid, grids(ilevel)%pw) 1576 END DO 1577 force(i) = ( energy(2) - energy(1) ) / (2.0_dp * Dx) 1578 END DO 1579 Err = 0.0_dp 1580 IF (ALL(force /= 0.0_dp)) THEN 1581 Err(1) = (debug_force(1)-force(1))/force(1)*100.0_dp 1582 Err(2) = (debug_force(2)-force(2))/force(2)*100.0_dp 1583 Err(3) = (debug_force(3)-force(3))/force(3)*100.0_dp 1584 END IF 1585 IF (iw>0) & 1586 WRITE(iw,100)Icount, debug_force(1), force(1), Err(1),& 1587 debug_force(2), force(2), Err(2),& 1588 debug_force(3), force(3), Err(3) 1589 CPASSERT(ABS(Err(1))<=MaxErr) 1590 CPASSERT(ABS(Err(2))<=MaxErr) 1591 CPASSERT(ABS(Err(3))<=MaxErr) 1592 1593 IF (ASSOCIATED(xdat)) THEN 1594 DEALLOCATE (xdat) 1595 ENDIF 1596 IF (ASSOCIATED(ydat)) THEN 1597 DEALLOCATE (ydat) 1598 ENDIF 1599 IF (ASSOCIATED(zdat)) THEN 1600 DEALLOCATE (zdat) 1601 ENDIF 1602 1603 CALL pw_pools_give_back_pws(aug_pools,grids) 1604 CALL timestop(handle) 1605100 FORMAT("Collocation : ",I5,2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ") 1606 END SUBROUTINE debug_integrate_gf_rspace_NoPBC 1607 1608! ************************************************************************************************** 1609!> \brief Debugs qmmm_forces_with_gaussian_LG.. It may helps too ... ;-] 1610!> \param pgfs ... 1611!> \param aug_pools ... 1612!> \param rho ... 1613!> \param mm_charges ... 1614!> \param mm_atom_index ... 1615!> \param mm_particles ... 1616!> \param num_mm_atoms ... 1617!> \param coarser_grid_level ... 1618!> \param per_potentials ... 1619!> \param debug_force ... 1620!> \param para_env ... 1621!> \param mm_cell ... 1622!> \param dOmmOqm ... 1623!> \param iw ... 1624!> \param par_scheme ... 1625!> \param qmmm_spherical_cutoff ... 1626!> \param shells ... 1627!> \par History 1628!> 08.2004 created [tlaino] 1629!> \author Teodoro Laino 1630! ************************************************************************************************** 1631 SUBROUTINE debug_qmmm_forces_with_gauss_LG(pgfs,aug_pools, rho, mm_charges, mm_atom_index,& 1632 mm_particles, num_mm_atoms, coarser_grid_level, per_potentials,& 1633 debug_force, para_env, mm_cell,dOmmOqm,iw,par_scheme,qmmm_spherical_cutoff, shells) 1634 1635 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 1636 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 1637 TYPE(pw_type), POINTER :: rho 1638 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 1639 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1640 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 1641 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level 1642 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials 1643 REAL(KIND=dp), DIMENSION(:, :) :: debug_force 1644 TYPE(cp_para_env_type), POINTER :: para_env 1645 TYPE(cell_type), POINTER :: mm_cell 1646 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 1647 INTEGER, INTENT(IN) :: iw, par_scheme 1648 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff 1649 LOGICAL :: shells 1650 1651 CHARACTER(len=*), PARAMETER :: routineN = 'debug_qmmm_forces_with_gauss_LG', & 1652 routineP = moduleN//':'//routineN 1653 1654 INTEGER :: handle, I, igrid, IndMM, J, K, ngrids 1655 REAL(KIND=dp) :: Coord_save 1656 REAL(KIND=dp), DIMENSION(2) :: energy 1657 REAL(KIND=dp), DIMENSION(3) :: Err 1658 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Num_Forces 1659 TYPE(pw_p_type), DIMENSION(:), POINTER :: grids 1660 1661 ALLOCATE(Num_Forces(3,num_mm_atoms)) 1662 NULLIFY(grids) 1663 CALL timeset(routineN,handle) 1664 ngrids = SIZE(aug_pools) 1665 CALL pw_pools_create_pws(aug_pools,grids,use_data=REALDATA3D,in_space=REALSPACE) 1666 DO igrid=1,ngrids 1667 CALL pw_zero(grids(igrid)%pw) 1668 END DO 1669 Atoms: DO I = 1, num_mm_atoms 1670 IndMM = mm_atom_index(I) 1671 Coords: DO J = 1, 3 1672 Coord_save = mm_particles(IndMM)%r(J) 1673 energy = 0.0_dp 1674 Diff: DO K = 1, 2 1675 mm_particles(IndMM)%r(J) = Coord_save + (-1)**K * Dx 1676 CALL pw_zero(grids(coarser_grid_level)%pw) 1677 1678 CALL qmmm_elec_with_gaussian_LG (pgfs=pgfs,& 1679 cgrid=grids(coarser_grid_level)%pw,& 1680 mm_charges=mm_charges,& 1681 mm_atom_index=mm_atom_index,& 1682 mm_particles=mm_particles,& 1683 para_env=para_env,& 1684 per_potentials=per_potentials,& 1685 mm_cell=mm_cell,& 1686 dOmmOqm=dOmmOqm,& 1687 par_scheme=par_scheme,& 1688 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 1689 shells=shells) 1690 1691 energy(K) = pw_integral_ab ( rho, grids(coarser_grid_level)%pw) 1692 END DO Diff 1693 IF (iw>0)& 1694 WRITE(iw,'(A,I6,A,I3,A,2F15.9)')& 1695 "DEBUG LR:: MM Atom = ",IndMM," Coord = ",J," Energies (+/-) :: ",energy(2), energy(1) 1696 Num_Forces(J,I) = ( energy(2) - energy(1) ) / (2.0_dp * Dx) 1697 mm_particles(IndMM)%r(J) = Coord_save 1698 END DO Coords 1699 END DO Atoms 1700 1701 DO I = 1, num_mm_atoms 1702 IndMM = mm_atom_index(I) 1703 Err = 0.0_dp 1704 IF (ALL(Num_Forces /= 0.0_dp)) THEN 1705 Err(1) = (debug_force(1,I)-Num_Forces(1,I))/Num_Forces(1,I)*100.0_dp 1706 Err(2) = (debug_force(2,I)-Num_Forces(2,I))/Num_Forces(2,I)*100.0_dp 1707 Err(3) = (debug_force(3,I)-Num_Forces(3,I))/Num_Forces(3,I)*100.0_dp 1708 END IF 1709 IF (iw>0) & 1710 WRITE(iw,100)IndMM,debug_force(1,I),Num_Forces(1,I),Err(1),& 1711 debug_force(2,I),Num_Forces(2,I),Err(2),& 1712 debug_force(3,I),Num_Forces(3,I),Err(3) 1713 CPASSERT(ABS(Err(1))<=MaxErr) 1714 CPASSERT(ABS(Err(2))<=MaxErr) 1715 CPASSERT(ABS(Err(3))<=MaxErr) 1716 END DO 1717 1718 DEALLOCATE(Num_Forces) 1719 CALL pw_pools_give_back_pws(aug_pools,grids) 1720 CALL timestop(handle) 1721100 FORMAT("MM Atom LR : ",I5,2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ") 1722 END SUBROUTINE debug_qmmm_forces_with_gauss_LG 1723 1724! ************************************************************************************************** 1725!> \brief Debugs qmmm_forces_with_gaussian_LR.. It may helps too ... ;-] 1726!> \param pgfs ... 1727!> \param aug_pools ... 1728!> \param rho ... 1729!> \param mm_charges ... 1730!> \param mm_atom_index ... 1731!> \param mm_particles ... 1732!> \param num_mm_atoms ... 1733!> \param coarser_grid_level ... 1734!> \param potentials ... 1735!> \param debug_force ... 1736!> \param para_env ... 1737!> \param mm_cell ... 1738!> \param dOmmOqm ... 1739!> \param iw ... 1740!> \param par_scheme ... 1741!> \param qmmm_spherical_cutoff ... 1742!> \param shells ... 1743!> \par History 1744!> 08.2004 created [tlaino] 1745!> \author Teodoro Laino 1746! ************************************************************************************************** 1747 SUBROUTINE debug_qmmm_forces_with_gauss_LR(pgfs,aug_pools, rho, mm_charges, mm_atom_index,& 1748 mm_particles, num_mm_atoms, coarser_grid_level, potentials,& 1749 debug_force, para_env, mm_cell,dOmmOqm,iw, par_scheme, qmmm_spherical_cutoff, shells) 1750 1751 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 1752 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools 1753 TYPE(pw_type), POINTER :: rho 1754 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 1755 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 1756 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles 1757 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level 1758 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: Potentials 1759 REAL(KIND=dp), DIMENSION(:, :) :: debug_force 1760 TYPE(cp_para_env_type), POINTER :: para_env 1761 TYPE(cell_type), POINTER :: mm_cell 1762 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 1763 INTEGER, INTENT(IN) :: iw, par_scheme 1764 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff 1765 LOGICAL :: shells 1766 1767 CHARACTER(len=*), PARAMETER :: routineN = 'debug_qmmm_forces_with_gauss_LR', & 1768 routineP = moduleN//':'//routineN 1769 1770 INTEGER :: handle, I, igrid, IndMM, J, K, ngrids 1771 REAL(KIND=dp) :: Coord_save 1772 REAL(KIND=dp), DIMENSION(2) :: energy 1773 REAL(KIND=dp), DIMENSION(3) :: Err 1774 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Num_Forces 1775 TYPE(pw_p_type), DIMENSION(:), POINTER :: grids 1776 1777 ALLOCATE(Num_Forces(3,num_mm_atoms)) 1778 NULLIFY(grids) 1779 CALL timeset(routineN,handle) 1780 ngrids = SIZE(aug_pools) 1781 CALL pw_pools_create_pws(aug_pools,grids,use_data=REALDATA3D,in_space=REALSPACE) 1782 DO igrid=1,ngrids 1783 CALL pw_zero(grids(igrid)%pw) 1784 END DO 1785 Atoms: DO I = 1, num_mm_atoms 1786 IndMM = mm_atom_index(I) 1787 Coords: DO J = 1, 3 1788 Coord_save = mm_particles(IndMM)%r(J) 1789 energy = 0.0_dp 1790 Diff: DO K = 1, 2 1791 mm_particles(IndMM)%r(J) = Coord_save + (-1)**K * Dx 1792 CALL pw_zero(grids(coarser_grid_level)%pw) 1793 1794 CALL qmmm_elec_with_gaussian_LR (pgfs=pgfs,& 1795 grid=grids(coarser_grid_level)%pw,& 1796 mm_charges=mm_charges,& 1797 mm_atom_index=mm_atom_index,& 1798 mm_particles=mm_particles,& 1799 para_env=para_env,& 1800 potentials=potentials,& 1801 mm_cell=mm_cell,& 1802 dOmmOqm=dOmmOqm,& 1803 par_scheme=par_scheme,& 1804 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 1805 shells=shells) 1806 1807 energy(K) = pw_integral_ab ( rho, grids(coarser_grid_level)%pw) 1808 END DO Diff 1809 IF (iw>0)& 1810 WRITE(iw,'(A,I6,A,I3,A,2F15.9)')& 1811 "DEBUG LR:: MM Atom = ",IndMM," Coord = ",J," Energies (+/-) :: ",energy(2), energy(1) 1812 Num_Forces(J,I) = ( energy(2) - energy(1) ) / (2.0_dp * Dx) 1813 mm_particles(IndMM)%r(J) = Coord_save 1814 END DO Coords 1815 END DO Atoms 1816 1817 DO I = 1, num_mm_atoms 1818 IndMM = mm_atom_index(I) 1819 Err = 0.0_dp 1820 IF (ALL(Num_Forces(:,I) /= 0.0_dp)) THEN 1821 Err(1) = (debug_force(1,I)-Num_Forces(1,I))/Num_Forces(1,I)*100.0_dp 1822 Err(2) = (debug_force(2,I)-Num_Forces(2,I))/Num_Forces(2,I)*100.0_dp 1823 Err(3) = (debug_force(3,I)-Num_Forces(3,I))/Num_Forces(3,I)*100.0_dp 1824 END IF 1825 IF (iw>0)& 1826 WRITE(iw,100)IndMM,debug_force(1,I),Num_Forces(1,I),Err(1),& 1827 debug_force(2,I),Num_Forces(2,I),Err(2),& 1828 debug_force(3,I),Num_Forces(3,I),Err(3) 1829 CPASSERT(ABS(Err(1))<=MaxErr) 1830 CPASSERT(ABS(Err(2))<=MaxErr) 1831 CPASSERT(ABS(Err(3))<=MaxErr) 1832 END DO 1833 1834 DEALLOCATE(Num_Forces) 1835 CALL pw_pools_give_back_pws(aug_pools,grids) 1836 CALL timestop(handle) 1837100 FORMAT("MM Atom LR : ",I5,2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ",2F15.9," ( ",F7.2," ) ") 1838 END SUBROUTINE debug_qmmm_forces_with_gauss_LR 1839 1840END MODULE qmmm_gpw_forces 1841