|
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 ! ***************************************************************************** 00012 MODULE qmmm_gpw_energy 00013 USE cell_types, ONLY: cell_type,& 00014 pbc 00015 USE cp_output_handling, ONLY: cp_p_file,& 00016 cp_print_key_finished_output,& 00017 cp_print_key_should_output,& 00018 cp_print_key_unit_nr 00019 USE cp_para_types, ONLY: cp_para_env_type 00020 USE cp_subsys_types, ONLY: cp_subsys_get,& 00021 cp_subsys_type 00022 USE cube_utils, ONLY: cube_info_type 00023 USE f77_blas 00024 USE input_constants, ONLY: do_par_atom,& 00025 do_qmmm_coulomb,& 00026 do_qmmm_gauss,& 00027 do_qmmm_none,& 00028 do_qmmm_swave,& 00029 spline3_nopbc_interp,& 00030 spline3_pbc_interp 00031 USE input_section_types, ONLY: section_get_ivals,& 00032 section_vals_get_subs_vals,& 00033 section_vals_type,& 00034 section_vals_val_get 00035 USE kinds, ONLY: dp 00036 USE message_passing, ONLY: mp_sum,& 00037 mp_sync 00038 USE mm_collocate_potential, ONLY: collocate_gf_rspace_NoPBC 00039 USE particle_list_types, ONLY: particle_list_type 00040 USE particle_types, ONLY: particle_type 00041 USE pw_env_types, ONLY: pw_env_get,& 00042 pw_env_type 00043 USE pw_methods, ONLY: pw_zero 00044 USE pw_pool_types, ONLY: pw_pool_p_type,& 00045 pw_pools_create_pws,& 00046 pw_pools_give_back_pws 00047 USE pw_spline_utils, ONLY: pw_prolongate_s3 00048 USE pw_types, ONLY: REALDATA3D,& 00049 REALSPACE,& 00050 pw_p_type,& 00051 pw_type 00052 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,& 00053 qmmm_gaussian_type 00054 USE qmmm_se_energy, ONLY: build_se_qmmm_matrix 00055 USE qmmm_types, ONLY: qmmm_env_qm_type,& 00056 qmmm_per_pot_p_type,& 00057 qmmm_per_pot_type,& 00058 qmmm_pot_p_type,& 00059 qmmm_pot_type 00060 USE qmmm_util, ONLY: spherical_cutoff_factor 00061 USE qs_environment_types, ONLY: get_qs_env,& 00062 qs_environment_type 00063 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type 00064 USE realspace_grid_cube, ONLY: pw_to_cube 00065 USE timings, ONLY: timeset,& 00066 timestop 00067 #include "cp_common_uses.h" 00068 00069 IMPLICIT NONE 00070 PRIVATE 00071 00072 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. 00073 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_energy' 00074 00075 PUBLIC :: qmmm_el_coupling 00076 PUBLIC :: qmmm_elec_with_gaussian, qmmm_elec_with_gaussian_low,& 00077 qmmm_elec_with_gaussian_LR, qmmm_elec_with_gaussian_LG 00078 !*** 00079 CONTAINS 00080 00081 ! ***************************************************************************** 00089 SUBROUTINE qmmm_el_coupling(qs_env,qmmm_env,mm_particles,mm_cell,error) 00090 TYPE(qs_environment_type), POINTER :: qs_env 00091 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 00092 TYPE(particle_type), DIMENSION(:), 00093 POINTER :: mm_particles 00094 TYPE(cell_type), POINTER :: mm_cell 00095 TYPE(cp_error_type), INTENT(inout) :: error 00096 00097 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_el_coupling', 00098 routineP = moduleN//':'//routineN 00099 00100 INTEGER :: handle, iw, iw2 00101 LOGICAL :: failure 00102 TYPE(cp_logger_type), POINTER :: logger 00103 TYPE(cp_para_env_type), POINTER :: para_env 00104 TYPE(cp_subsys_type), POINTER :: subsys 00105 TYPE(particle_list_type), POINTER :: particles 00106 TYPE(pw_env_type), POINTER :: pw_env 00107 TYPE(pw_pool_p_type), DIMENSION(:), 00108 POINTER :: pw_pools 00109 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc 00110 TYPE(section_vals_type), POINTER :: input_section, 00111 interp_section, print_section 00112 00113 CALL timeset(routineN,handle) 00114 failure=.FALSE. 00115 logger => cp_error_get_logger(error) 00116 IF (.NOT.failure) THEN 00117 NULLIFY(ks_qmmm_env_loc, pw_pools, pw_env,input_section) 00118 CALL get_qs_env(qs_env=qs_env,& 00119 pw_env=pw_env,& 00120 para_env=para_env,& 00121 input=input_section,& 00122 ks_qmmm_env=ks_qmmm_env_loc,& 00123 subsys=subsys,& 00124 error=error) 00125 CALL cp_subsys_get(subsys,particles=particles,error=error) 00126 00127 00128 CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools, error=error) 00129 print_section => section_vals_get_subs_vals(input_section,"QMMM%PRINT",error=error) 00130 iw = cp_print_key_unit_nr(logger,print_section,"PROGRAM_RUN_INFO",& 00131 extension=".qmmmLog",error=error) 00132 IF (iw>0) & 00133 WRITE(iw,'(T2,"QMMM|",1X,A)')"Information on the QM/MM Electrostatic Potential:" 00134 ! 00135 ! Initializing vectors: 00136 ! Zeroing v_qmmm_rspace 00137 CALL pw_zero(ks_qmmm_env_loc%v_qmmm_rspace%pw,error=error) 00138 IF (qs_env%dft_control%qs_control%semi_empirical) THEN 00139 ! SEMIEMPIRICAL 00140 SELECT CASE(qmmm_env%qmmm_coupl_type) 00141 CASE(do_qmmm_coulomb,do_qmmm_none) 00142 CALL build_se_qmmm_matrix(qs_env,qmmm_env,mm_particles,mm_cell,qs_env%para_env,error) 00143 IF( qmmm_env%qmmm_coupl_type==do_qmmm_none) THEN 00144 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')& 00145 "No QM/MM Electrostatic coupling. Just Mechanical Coupling!" 00146 END IF 00147 CASE (do_qmmm_gauss,do_qmmm_swave) 00148 CALL cp_unimplemented_error(fromWhere=routineP, & 00149 message="GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.", & 00150 error=error, error_level=cp_failure_level) 00151 CASE DEFAULT 00152 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')"Unknown Coupling..." 00153 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00154 END SELECT 00155 ELSEIF (qs_env%dft_control%qs_control%dftb) THEN 00156 ! DFTB 00157 CALL cp_unimplemented_error(fromWhere=routineP, & 00158 message="QM/MM electrostatic coupling not yet implemented for DFTB.", & 00159 error=error, error_level=cp_failure_level) 00160 ELSE 00161 ! QS 00162 SELECT CASE(qmmm_env%qmmm_coupl_type) 00163 CASE(do_qmmm_coulomb) 00164 CALL cp_unimplemented_error(fromWhere=routineP, & 00165 message="Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.", & 00166 error=error, error_level=cp_failure_level) 00167 CASE(do_qmmm_gauss,do_qmmm_swave) 00168 IF (iw>0) & 00169 WRITE(iw,'(T2,"QMMM|",1X,A)')& 00170 "QM/MM Coupling computed collocating the Gaussian Potential Functions." 00171 interp_section => section_vals_get_subs_vals(input_section,& 00172 "QMMM%INTERPOLATOR",error=error) 00173 CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env,& 00174 v_qmmm=ks_qmmm_env_loc%v_qmmm_rspace,& 00175 mm_particles=mm_particles,& 00176 aug_pools=qmmm_env%aug_pools,& 00177 para_env=para_env,& 00178 eps_mm_rspace=qmmm_env%eps_mm_rspace,& 00179 cube_info=ks_qmmm_env_loc%cube_info,& 00180 pw_pools=pw_pools,& 00181 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid,& 00182 coarser_grid=qmmm_env%gridlevel_info%coarser_grid,& 00183 interp_section=interp_section,& 00184 mm_cell=mm_cell,& 00185 error=error) 00186 CASE(do_qmmm_none) 00187 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')& 00188 "No QM/MM Electrostatic coupling. Just Mechanical Coupling!" 00189 CASE DEFAULT 00190 IF (iw>0) WRITE(iw,'(T2,"QMMM|",1X,A)')"Unknown Coupling..." 00191 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00192 END SELECT 00193 ! Dump info on the electrostatic potential if requested 00194 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_section,& 00195 "POTENTIAL",error=error),cp_p_file)) THEN 00196 iw2 = cp_print_key_unit_nr(logger,print_section,"POTENTIAL",& 00197 extension=".qmmmLog",error=error) 00198 CALL pw_to_cube(ks_qmmm_env_loc%v_qmmm_rspace%pw,iw2,& 00199 particles=particles,& 00200 stride=section_get_ivals(print_section,"POTENTIAL%STRIDE",error),& 00201 title="QM/MM: MM ELECTROSTATIC POTENTIAL ", error=error) 00202 CALL cp_print_key_finished_output(iw2,logger,print_section,& 00203 "POTENTIAL", error=error) 00204 END IF 00205 END IF 00206 CALL cp_print_key_finished_output(iw,logger,print_section,& 00207 "PROGRAM_RUN_INFO", error=error) 00208 CALL timestop(handle) 00209 END IF 00210 END SUBROUTINE qmmm_el_coupling 00211 00212 ! ***************************************************************************** 00221 SUBROUTINE qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles,& 00222 aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools,& 00223 auxbas_grid, coarser_grid, interp_section, mm_cell, error) 00224 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 00225 TYPE(pw_p_type), INTENT(INOUT) :: v_qmmm 00226 TYPE(particle_type), DIMENSION(:), 00227 POINTER :: mm_particles 00228 TYPE(pw_pool_p_type), DIMENSION(:), 00229 POINTER :: aug_pools 00230 TYPE(cube_info_type), DIMENSION(:), 00231 POINTER :: cube_info 00232 TYPE(cp_para_env_type), POINTER :: para_env 00233 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace 00234 TYPE(pw_pool_p_type), DIMENSION(:), 00235 POINTER :: pw_pools 00236 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid 00237 TYPE(section_vals_type), POINTER :: interp_section 00238 TYPE(cell_type), POINTER :: mm_cell 00239 TYPE(cp_error_type), INTENT(inout) :: error 00240 00241 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian', 00242 routineP = moduleN//':'//routineN 00243 00244 INTEGER :: handle, handle2, igrid, 00245 ilevel, kind_interp, lb(3), 00246 ngrids, ub(3) 00247 LOGICAL :: failure 00248 TYPE(pw_p_type), DIMENSION(:), POINTER :: grids 00249 00250 failure=.FALSE. 00251 CPPrecondition(ASSOCIATED(mm_particles),cp_failure_level,routineP,error,failure) 00252 CPPrecondition(ASSOCIATED(qmmm_env%mm_atom_chrg),cp_failure_level,routineP,error,failure) 00253 CPPrecondition(ASSOCIATED(qmmm_env%mm_atom_index),cp_failure_level,routineP,error,failure) 00254 CPPrecondition(ASSOCIATED(aug_pools),cp_failure_level,routineP,error,failure) 00255 CPPrecondition(ASSOCIATED(pw_pools),cp_failure_level,routineP,error,failure) 00256 !Statements 00257 CALL timeset(routineN,handle) 00258 ngrids=SIZE(pw_pools) 00259 CALL pw_pools_create_pws(aug_pools,grids,use_data=REALDATA3D,in_space=REALSPACE,error=error) 00260 DO igrid=1,ngrids 00261 CALL pw_zero(grids(igrid)%pw,error=error) 00262 END DO 00263 00264 CALL qmmm_elec_with_gaussian_low( grids, mm_particles,& 00265 qmmm_env%mm_atom_chrg, qmmm_env%mm_el_pot_radius, qmmm_env%mm_atom_index, & 00266 qmmm_env%num_mm_atoms, cube_info, para_env, eps_mm_rspace, qmmm_env%pgfs, & 00267 auxbas_grid, coarser_grid, aug_pools, qmmm_env%potentials, & 00268 mm_cell=mm_cell, dOmmOqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, & 00269 per_potentials=qmmm_env%per_potentials, par_scheme=qmmm_env%par_scheme, & 00270 qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, error=error) 00271 00272 IF (qmmm_env%move_mm_charges.OR.qmmm_env%add_mm_charges) THEN 00273 CALL qmmm_elec_with_gaussian_low( grids, qmmm_env%added_charges%added_particles, & 00274 qmmm_env%added_charges%mm_atom_chrg, qmmm_env%added_charges%mm_el_pot_radius, & 00275 qmmm_env%added_charges%mm_atom_index, qmmm_env%added_charges%num_mm_atoms, & 00276 cube_info, para_env, eps_mm_rspace, qmmm_env%added_charges%pgfs,auxbas_grid, & 00277 coarser_grid, aug_pools, qmmm_env%added_charges%potentials, & 00278 mm_cell=mm_cell, dOmmOqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, & 00279 per_potentials=qmmm_env%per_potentials, par_scheme=qmmm_env%par_scheme, & 00280 qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, error=error) 00281 END IF 00282 ! Sumup all contributions according the parallelization scheme 00283 IF (qmmm_env%par_scheme==do_par_atom) THEN 00284 DO ilevel=1,SIZE(grids) 00285 CALL mp_sum(grids(ilevel)%pw%cr3d,para_env%group) 00286 END DO 00287 END IF 00288 ! RealSpace Interpolation 00289 CALL section_vals_val_get(interp_section,"kind", i_val=kind_interp, error=error) 00290 SELECT CASE(kind_interp) 00291 CASE(spline3_nopbc_interp, spline3_pbc_interp) 00292 ! Spline Iterpolator 00293 CALL mp_sync(para_env%group) 00294 CALL timeset(TRIM(routineN)//":spline3Int",handle2) 00295 DO Ilevel = coarser_grid, auxbas_grid+1, -1 00296 CALL pw_prolongate_s3(grids(Ilevel )%pw,& 00297 grids(Ilevel-1)%pw,& 00298 aug_pools(Ilevel)%pool,& 00299 param_section=interp_section,& 00300 error=error) 00301 END DO 00302 CALL timestop(handle2) 00303 CASE DEFAULT 00304 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00305 END SELECT 00306 lb = v_qmmm%pw%pw_grid%bounds_local(1,:) 00307 ub = v_qmmm%pw%pw_grid%bounds_local(2,:) 00308 00309 v_qmmm%pw%cr3d = grids(auxbas_grid)%pw%cr3d (lb(1):ub(1),& 00310 lb(2):ub(2),& 00311 lb(3):ub(3) ) 00312 00313 CALL pw_pools_give_back_pws(aug_pools,grids,error=error) 00314 00315 CALL timestop(handle) 00316 END SUBROUTINE qmmm_elec_with_gaussian 00317 00318 ! ***************************************************************************** 00327 SUBROUTINE qmmm_elec_with_gaussian_low( tmp_grid, mm_particles, mm_charges,& 00328 mm_el_pot_radius, mm_atom_index, num_mm_atoms, cube_info, para_env, & 00329 eps_mm_rspace, pgfs, auxbas_grid, coarser_grid, aug_pools, & 00330 potentials, mm_cell, dOmmOqm, periodic, per_potentials, par_scheme, & 00331 qmmm_spherical_cutoff, error) 00332 TYPE(pw_p_type), DIMENSION(:), POINTER :: tmp_grid 00333 TYPE(particle_type), DIMENSION(:), 00334 POINTER :: mm_particles 00335 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges, mm_el_pot_radius 00336 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 00337 INTEGER, INTENT(IN) :: num_mm_atoms 00338 TYPE(cube_info_type), DIMENSION(:), 00339 POINTER :: cube_info 00340 TYPE(cp_para_env_type), POINTER :: para_env 00341 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace 00342 TYPE(qmmm_gaussian_p_type), 00343 DIMENSION(:), POINTER :: pgfs 00344 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid 00345 TYPE(pw_pool_p_type), DIMENSION(:), 00346 POINTER :: aug_pools 00347 TYPE(qmmm_pot_p_type), DIMENSION(:), 00348 POINTER :: potentials 00349 TYPE(cell_type), POINTER :: mm_cell 00350 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 00351 LOGICAL, INTENT(IN) :: periodic 00352 TYPE(qmmm_per_pot_p_type), 00353 DIMENSION(:), POINTER :: per_potentials 00354 INTEGER, INTENT(IN) :: par_scheme 00355 REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2) 00356 TYPE(cp_error_type), INTENT(inout) :: error 00357 00358 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian_low', 00359 routineNb = 'qmmm_elec_gaussian_low', routineP = moduleN//':'//routineN 00360 00361 INTEGER :: handle, handle2, IGauss, 00362 ilevel, Imm, IndMM, IRadTyp, 00363 LIndMM, myind, n_rep_real(3), 00364 stat 00365 INTEGER, DIMENSION(2, 3) :: bo2 00366 LOGICAL :: failure 00367 REAL(KIND=dp) :: alpha, height, 00368 sph_chrg_factor, W 00369 REAL(KIND=dp), DIMENSION(3) :: ra 00370 REAL(KIND=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat 00371 TYPE(qmmm_gaussian_type), POINTER :: pgf 00372 TYPE(qmmm_per_pot_type), POINTER :: per_pot 00373 TYPE(qmmm_pot_type), POINTER :: pot 00374 00375 NULLIFY(pgf, pot, per_pot, xdat, ydat, zdat) 00376 CALL timeset(routineN,handle) 00377 CALL timeset(routineNb//"_G",handle2) 00378 bo2 = tmp_grid(auxbas_grid)%pw%pw_grid%bounds 00379 ALLOCATE (xdat(bo2(1,1):bo2(2,1)), STAT=stat) 00380 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00381 ALLOCATE (ydat(bo2(1,2):bo2(2,2)), STAT=stat) 00382 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00383 ALLOCATE (zdat(bo2(1,3):bo2(2,3)), STAT=stat) 00384 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00385 IF (par_scheme==do_par_atom) myind = 0 00386 Radius: DO IRadTyp = 1, SIZE(pgfs) 00387 pgf => pgfs(IRadTyp)%pgf 00388 pot => potentials(IRadTyp)%pot 00389 n_rep_real = 0 00390 IF (periodic) THEN 00391 per_pot => per_potentials(IRadTyp)%pot 00392 n_rep_real = per_pot%n_rep_real 00393 END IF 00394 Gaussian: DO IGauss = 1, pgf%Number_of_Gaussians 00395 alpha = 1.0_dp / pgf%Gk(IGauss) 00396 alpha = alpha * alpha 00397 height = pgf%Ak(IGauss) 00398 ilevel = pgf%grid_level(IGauss) 00399 Atoms: DO Imm = 1, SIZE(pot%mm_atom_index) 00400 IF (par_scheme==do_par_atom) THEN 00401 myind = myind + 1 00402 IF (MOD(myind,para_env%num_pe)/=para_env%mepos) CYCLE Atoms 00403 END IF 00404 LIndMM = pot%mm_atom_index(Imm) 00405 IndMM = mm_atom_index(LIndMM) 00406 ra(:) = pbc(mm_particles(IndMM)%r-dOmmOqm, mm_cell)+dOmmOqm 00407 W = mm_charges(LIndMM) * height 00408 ! Possible Spherical Cutoff 00409 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 00410 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor, error) 00411 W = W * sph_chrg_factor 00412 END IF 00413 IF (ABS(W)<= EPSILON(0.0_dp)) CYCLE Atoms 00414 CALL collocate_gf_rspace_NoPBC(zetp=alpha,& 00415 rp=ra,& 00416 scale=-1.0_dp,& 00417 W=W,& 00418 pwgrid=tmp_grid(ilevel)%pw,& 00419 cube_info=cube_info(ilevel),& 00420 eps_mm_rspace=eps_mm_rspace,& 00421 xdat=xdat,& 00422 ydat=ydat,& 00423 zdat=zdat,& 00424 bo2=bo2,& 00425 n_rep_real=n_rep_real,& 00426 mm_cell=mm_cell) 00427 END DO Atoms 00428 END DO Gaussian 00429 END DO Radius 00430 IF (ASSOCIATED(xdat)) THEN 00431 DEALLOCATE (xdat, STAT=stat) 00432 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00433 ENDIF 00434 IF (ASSOCIATED(ydat)) THEN 00435 DEALLOCATE (ydat, STAT=stat) 00436 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00437 ENDIF 00438 IF (ASSOCIATED(zdat)) THEN 00439 DEALLOCATE (zdat, STAT=stat) 00440 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00441 ENDIF 00442 CALL timestop(handle2) 00443 CALL timeset(routineNb//"_R",handle2) 00444 IF (periodic) THEN 00445 ! Long Range Part of the QM/MM Potential with Gaussians With Periodic Boundary Conditions 00446 CALL qmmm_elec_with_gaussian_LG (pgfs=pgfs,& 00447 cgrid=tmp_grid(coarser_grid)%pw,& 00448 mm_charges=mm_charges,& 00449 mm_atom_index=mm_atom_index,& 00450 mm_particles=mm_particles,& 00451 para_env=para_env,& 00452 coarser_grid_level=coarser_grid,& 00453 per_potentials=per_potentials,& 00454 mm_cell=mm_cell,& 00455 dOmmOqm=dOmmOqm,& 00456 par_scheme=par_scheme,& 00457 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 00458 error=error) 00459 ELSE 00460 ! Long Range Part of the QM/MM Potential with Gaussians 00461 CALL qmmm_elec_with_gaussian_LR (pgfs=pgfs,& 00462 grid=tmp_grid(coarser_grid)%pw,& 00463 mm_charges=mm_charges,& 00464 mm_atom_index=mm_atom_index,& 00465 mm_particles=mm_particles,& 00466 mm_el_pot_radius=mm_el_pot_radius,& 00467 para_env=para_env,& 00468 coarser_grid_level=coarser_grid,& 00469 potentials=potentials,& 00470 mm_cell=mm_cell,& 00471 dOmmOqm=dOmmOqm,& 00472 par_scheme=par_scheme,& 00473 qmmm_spherical_cutoff=qmmm_spherical_cutoff,& 00474 error=error) 00475 END IF 00476 CALL timestop(handle2) 00477 CALL timestop(handle) 00478 00479 END SUBROUTINE qmmm_elec_with_gaussian_low 00480 00481 ! ***************************************************************************** 00495 SUBROUTINE qmmm_elec_with_gaussian_LG(pgfs, cgrid, mm_charges, mm_atom_index,& 00496 mm_particles, para_env, coarser_grid_level, per_potentials,& 00497 mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, error) 00498 TYPE(qmmm_gaussian_p_type), 00499 DIMENSION(:), POINTER :: pgfs 00500 TYPE(pw_type), POINTER :: cgrid 00501 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 00502 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 00503 TYPE(particle_type), DIMENSION(:), 00504 POINTER :: mm_particles 00505 TYPE(cp_para_env_type), POINTER :: para_env 00506 INTEGER, INTENT(IN) :: coarser_grid_level 00507 TYPE(qmmm_per_pot_p_type), 00508 DIMENSION(:), POINTER :: per_potentials 00509 TYPE(cell_type), POINTER :: mm_cell 00510 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 00511 INTEGER, INTENT(IN) :: par_scheme 00512 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff 00513 TYPE(cp_error_type), INTENT(inout) :: error 00514 00515 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian_LG', 00516 routineP = moduleN//':'//routineN 00517 00518 INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, ij3, ij4, ik1, ik2, 00519 ik3, ik4, Imm, IndMM, IRadTyp, ivec(3), j, k, LIndMM, my_j, my_k, 00520 myind, npts(3) 00521 INTEGER, DIMENSION(2, 3) :: bo, gbo 00522 LOGICAL :: failure 00523 REAL(KIND=dp) :: a1, a2, a3, abc_X(4,4), abc_X_Y(4), b1, b2, b3, c1, c2, 00524 c3, d1, d2, d3, dr1, dr1c, dr2, dr2c, dr3, dr3c, e1, e2, e3, f1, f2, 00525 f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, 00526 rt1, rt2, rt3, rv1, rv2, rv3, s1, s2, s3, s4, sph_chrg_factor, t1, t2, 00527 t3, t4, u1, u2, u3, v1, v2, v3, v4, val, xd1, xd2, xd3, xs1, xs2, xs3 00528 REAL(KIND=dp), DIMENSION(3) :: ra, vec 00529 REAL(KIND=dp), DIMENSION(:, :, :), 00530 POINTER :: grid, grid2 00531 TYPE(pw_type), POINTER :: pw 00532 TYPE(qmmm_per_pot_type), POINTER :: per_pot 00533 00534 failure = .FALSE. 00535 CALL timeset(routineN,handle) 00536 NULLIFY(grid, pw) 00537 dr1c = cgrid%pw_grid%dr(1) 00538 dr2c = cgrid%pw_grid%dr(2) 00539 dr3c = cgrid%pw_grid%dr(3) 00540 gbo = cgrid%pw_grid%bounds 00541 bo = cgrid%pw_grid%bounds_local 00542 grid2=>cgrid%cr3d 00543 IF (par_scheme==do_par_atom) myind = 0 00544 Radius: DO IRadTyp = 1, SIZE(pgfs) 00545 per_pot => per_potentials(IRadTyp)%pot 00546 pw => per_pot%TabLR 00547 npts = pw%pw_grid%npts 00548 dr1 = pw%pw_grid%dr(1) 00549 dr2 = pw%pw_grid%dr(2) 00550 dr3 = pw%pw_grid%dr(3) 00551 grid => pw%cr3d(:,:,:) 00552 Atoms: DO Imm = 1, SIZE(per_pot%mm_atom_index) 00553 IF (par_scheme==do_par_atom) THEN 00554 myind = myind + 1 00555 IF (MOD(myind,para_env%num_pe)/=para_env%mepos) CYCLE 00556 END IF 00557 LIndMM = per_pot%mm_atom_index(Imm) 00558 IndMM = mm_atom_index(LIndMM) 00559 ra(:) = pbc(mm_particles(IndMM)%r-dOmmOqm,mm_cell)+dOmmOqm 00560 qt = mm_charges(LIndMM) 00561 ! Possible Spherical Cutoff 00562 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 00563 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor, error) 00564 qt = qt * sph_chrg_factor 00565 END IF 00566 IF (ABS(qt)<= EPSILON(0.0_dp)) CYCLE Atoms 00567 rt1 = ra(1) 00568 rt2 = ra(2) 00569 rt3 = ra(3) 00570 LoopOnGrid: DO k = bo(1,3), bo(2,3) 00571 my_k=k-gbo(1,3) 00572 xs3 = REAL(my_k,dp)*dr3c 00573 my_j=bo(1,2)-gbo(1,2) 00574 xs2 = REAL(my_j,dp)*dr2c 00575 rv3 = rt3 - xs3 00576 vec(3) = rv3 00577 ivec(3) = FLOOR(vec(3)/pw%pw_grid%dr(3)) 00578 xd3 = (vec(3)/dr3)-REAL(ivec(3),kind=dp) 00579 ik1 = MODULO(ivec(3)-1,npts(3))+1 00580 ik2 = MODULO(ivec(3) ,npts(3))+1 00581 ik3 = MODULO(ivec(3)+1,npts(3))+1 00582 ik4 = MODULO(ivec(3)+2,npts(3))+1 00583 p1 = 3.0_dp + xd3 00584 p2 = p1*p1 00585 p3 = p2*p1 00586 q1 = 2.0_dp + xd3 00587 q2 = q1*q1 00588 q3 = q2*q1 00589 r1 = 1.0_dp + xd3 00590 r2 = r1*r1 00591 r3 = r2*r1 00592 u1 = xd3 00593 u2 = u1*u1 00594 u3 = u2*u1 00595 v1 = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3) 00596 v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3 00597 v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3 00598 v4 = 1.0_dp/6.0_dp*u3 00599 DO j = bo(1,2), bo(2,2) 00600 xs1 = (bo(1,1)-gbo(1,1))*dr1c 00601 rv2 = rt2 - xs2 00602 vec(2) = rv2 00603 ivec(2) = FLOOR(vec(2)/pw%pw_grid%dr(2)) 00604 xd2 = (vec(2)/dr2)-REAL(ivec(2),kind=dp) 00605 ij1 = MODULO(ivec(2)-1,npts(2))+1 00606 ij2 = MODULO(ivec(2) ,npts(2))+1 00607 ij3 = MODULO(ivec(2)+1,npts(2))+1 00608 ij4 = MODULO(ivec(2)+2,npts(2))+1 00609 e1 = 3.0_dp + xd2 00610 e2 = e1*e1 00611 e3 = e2*e1 00612 f1 = 2.0_dp + xd2 00613 f2 = f1*f1 00614 f3 = f2*f1 00615 g1 = 1.0_dp + xd2 00616 g2 = g1*g1 00617 g3 = g2*g1 00618 h1 = xd2 00619 h2 = h1*h1 00620 h3 = h2*h1 00621 s1 = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3) 00622 s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3 00623 s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3 00624 s4 = 1.0_dp/6.0_dp*h3 00625 DO i = bo(1,1), bo(2,1) 00626 rv1 = rt1 - xs1 00627 vec(1) = rv1 00628 ivec(1) = FLOOR(vec(1)/pw%pw_grid%dr(1)) 00629 xd1 = (vec(1)/dr1)-REAL(ivec(1),kind=dp) 00630 ii1 = MODULO(ivec(1)-1,npts(1))+1 00631 ii2 = MODULO(ivec(1) ,npts(1))+1 00632 ii3 = MODULO(ivec(1)+1,npts(1))+1 00633 ii4 = MODULO(ivec(1)+2,npts(1))+1 00634 ! 00635 ! Spline Interpolation 00636 ! 00637 00638 a1 = 3.0_dp + xd1 00639 a2 = a1*a1 00640 a3 = a2*a1 00641 b1 = 2.0_dp + xd1 00642 b2 = b1*b1 00643 b3 = b2*b1 00644 c1 = 1.0_dp + xd1 00645 c2 = c1*c1 00646 c3 = c2*c1 00647 d1 = xd1 00648 d2 = d1*d1 00649 d3 = d2*d1 00650 t1 = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3) 00651 t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3 00652 t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3 00653 t4 = 1.0_dp/6.0_dp*d3 00654 00655 abc_X(1,1) = grid(ii1,ij1,ik1)*v1 + grid(ii1,ij1,ik2)*v2 + grid(ii1,ij1,ik3)*v3 + grid(ii1,ij1,ik4)*v4 00656 abc_X(1,2) = grid(ii1,ij2,ik1)*v1 + grid(ii1,ij2,ik2)*v2 + grid(ii1,ij2,ik3)*v3 + grid(ii1,ij2,ik4)*v4 00657 abc_X(1,3) = grid(ii1,ij3,ik1)*v1 + grid(ii1,ij3,ik2)*v2 + grid(ii1,ij3,ik3)*v3 + grid(ii1,ij3,ik4)*v4 00658 abc_X(1,4) = grid(ii1,ij4,ik1)*v1 + grid(ii1,ij4,ik2)*v2 + grid(ii1,ij4,ik3)*v3 + grid(ii1,ij4,ik4)*v4 00659 abc_X(2,1) = grid(ii2,ij1,ik1)*v1 + grid(ii2,ij1,ik2)*v2 + grid(ii2,ij1,ik3)*v3 + grid(ii2,ij1,ik4)*v4 00660 abc_X(2,2) = grid(ii2,ij2,ik1)*v1 + grid(ii2,ij2,ik2)*v2 + grid(ii2,ij2,ik3)*v3 + grid(ii2,ij2,ik4)*v4 00661 abc_X(2,3) = grid(ii2,ij3,ik1)*v1 + grid(ii2,ij3,ik2)*v2 + grid(ii2,ij3,ik3)*v3 + grid(ii2,ij3,ik4)*v4 00662 abc_X(2,4) = grid(ii2,ij4,ik1)*v1 + grid(ii2,ij4,ik2)*v2 + grid(ii2,ij4,ik3)*v3 + grid(ii2,ij4,ik4)*v4 00663 abc_X(3,1) = grid(ii3,ij1,ik1)*v1 + grid(ii3,ij1,ik2)*v2 + grid(ii3,ij1,ik3)*v3 + grid(ii3,ij1,ik4)*v4 00664 abc_X(3,2) = grid(ii3,ij2,ik1)*v1 + grid(ii3,ij2,ik2)*v2 + grid(ii3,ij2,ik3)*v3 + grid(ii3,ij2,ik4)*v4 00665 abc_X(3,3) = grid(ii3,ij3,ik1)*v1 + grid(ii3,ij3,ik2)*v2 + grid(ii3,ij3,ik3)*v3 + grid(ii3,ij3,ik4)*v4 00666 abc_X(3,4) = grid(ii3,ij4,ik1)*v1 + grid(ii3,ij4,ik2)*v2 + grid(ii3,ij4,ik3)*v3 + grid(ii3,ij4,ik4)*v4 00667 abc_X(4,1) = grid(ii4,ij1,ik1)*v1 + grid(ii4,ij1,ik2)*v2 + grid(ii4,ij1,ik3)*v3 + grid(ii4,ij1,ik4)*v4 00668 abc_X(4,2) = grid(ii4,ij2,ik1)*v1 + grid(ii4,ij2,ik2)*v2 + grid(ii4,ij2,ik3)*v3 + grid(ii4,ij2,ik4)*v4 00669 abc_X(4,3) = grid(ii4,ij3,ik1)*v1 + grid(ii4,ij3,ik2)*v2 + grid(ii4,ij3,ik3)*v3 + grid(ii4,ij3,ik4)*v4 00670 abc_X(4,4) = grid(ii4,ij4,ik1)*v1 + grid(ii4,ij4,ik2)*v2 + grid(ii4,ij4,ik3)*v3 + grid(ii4,ij4,ik4)*v4 00671 00672 abc_X_Y(1) = abc_X(1,1)*t1 + abc_X(2,1)*t2 + abc_X(3,1)*t3 + abc_X(4,1)*t4 00673 abc_X_Y(2) = abc_X(1,2)*t1 + abc_X(2,2)*t2 + abc_X(3,2)*t3 + abc_X(4,2)*t4 00674 abc_X_Y(3) = abc_X(1,3)*t1 + abc_X(2,3)*t2 + abc_X(3,3)*t3 + abc_X(4,3)*t4 00675 abc_X_Y(4) = abc_X(1,4)*t1 + abc_X(2,4)*t2 + abc_X(3,4)*t3 + abc_X(4,4)*t4 00676 00677 val = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4 00678 00679 grid2(i,j,k) = grid2(i,j,k) - val * qt 00680 xs1 = xs1 + dr1c 00681 END DO 00682 xs2 = xs2 + dr2c 00683 END DO 00684 END DO LoopOnGrid 00685 END DO Atoms 00686 END DO Radius 00687 CALL timestop(handle) 00688 END SUBROUTINE qmmm_elec_with_gaussian_LG 00689 00690 ! ***************************************************************************** 00700 SUBROUTINE qmmm_elec_with_gaussian_LR(pgfs, grid, mm_charges, mm_atom_index,& 00701 mm_particles, mm_el_pot_radius, para_env,coarser_grid_level, potentials,& 00702 mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, error) 00703 TYPE(qmmm_gaussian_p_type), 00704 DIMENSION(:), POINTER :: pgfs 00705 TYPE(pw_type), POINTER :: grid 00706 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 00707 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 00708 TYPE(particle_type), DIMENSION(:), 00709 POINTER :: mm_particles 00710 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_el_pot_radius 00711 TYPE(cp_para_env_type), POINTER :: para_env 00712 INTEGER, INTENT(IN) :: coarser_grid_level 00713 TYPE(qmmm_pot_p_type), DIMENSION(:), 00714 POINTER :: potentials 00715 TYPE(cell_type), POINTER :: mm_cell 00716 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dOmmOqm 00717 INTEGER, INTENT(IN) :: par_scheme 00718 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff 00719 TYPE(cp_error_type), INTENT(inout) :: error 00720 00721 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian_LR', 00722 routineP = moduleN//':'//routineN 00723 00724 INTEGER :: handle, i, Imm, IndMM, 00725 IRadTyp, ix, j, k, LIndMM, 00726 my_j, my_k, myind, n1, n2, n3 00727 INTEGER, DIMENSION(2, 3) :: bo, gbo 00728 LOGICAL :: failure 00729 REAL(KIND=dp) :: dr1, dr2, dr3, dx, qt, r, r2, 00730 rt1, rt2, rt3, rv1, rv2, rv3, 00731 rx, rx2, rx3, 00732 sph_chrg_factor, Term, xs1, 00733 xs2, xs3 00734 REAL(KIND=dp), DIMENSION(3) :: ra 00735 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pot0_2 00736 REAL(KIND=dp), DIMENSION(:, :, :), 00737 POINTER :: grid2 00738 TYPE(qmmm_pot_type), POINTER :: pot 00739 00740 CALL timeset(routineN,handle) 00741 failure = .FALSE. 00742 n1 = grid%pw_grid%npts(1) 00743 n2 = grid%pw_grid%npts(2) 00744 n3 = grid%pw_grid%npts(3) 00745 dr1 = grid%pw_grid%dr(1) 00746 dr2 = grid%pw_grid%dr(2) 00747 dr3 = grid%pw_grid%dr(3) 00748 gbo = grid%pw_grid%bounds 00749 bo = grid%pw_grid%bounds_local 00750 grid2=>grid%cr3d 00751 IF (par_scheme==do_par_atom) myind=0 00752 Radius: DO IRadTyp = 1, SIZE(pgfs) 00753 pot => potentials(IRadTyp)%pot 00754 dx = Pot%dx 00755 pot0_2 => Pot%pot0_2 00756 Atoms: DO Imm = 1, SIZE(pot%mm_atom_index) 00757 IF (par_scheme==do_par_atom) THEN 00758 myind = myind + 1 00759 IF (MOD(myind,para_env%num_pe)/=para_env%mepos) CYCLE 00760 END IF 00761 LIndMM = pot%mm_atom_index(Imm) 00762 IndMM = mm_atom_index(LIndMM) 00763 ra(:) = pbc(mm_particles(IndMM)%r-dOmmOqm,mm_cell)+dOmmOqm 00764 qt = mm_charges(LIndMM) 00765 ! Possible Spherical Cutoff 00766 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 00767 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor, error) 00768 qt = qt * sph_chrg_factor 00769 END IF 00770 IF (ABS(qt)<= EPSILON(0.0_dp)) CYCLE Atoms 00771 rt1 = ra(1) 00772 rt2 = ra(2) 00773 rt3 = ra(3) 00774 LoopOnGrid: DO k = bo(1,3), bo(2,3) 00775 my_k=k-gbo(1,3) 00776 xs3 = REAL(my_k,dp)*dr3 00777 my_j=bo(1,2)-gbo(1,2) 00778 xs2 = REAL(my_j,dp)*dr2 00779 rv3 = rt3 - xs3 00780 DO j = bo(1,2), bo(2,2) 00781 xs1 = (bo(1,1)-gbo(1,1))*dr1 00782 rv2 = rt2 - xs2 00783 DO i = bo(1,1), bo(2,1) 00784 rv1 = rt1 - xs1 00785 r2 = rv1*rv1 + rv2*rv2 + rv3*rv3 00786 r = SQRT(r2) 00787 ix = FLOOR(r/dx)+1 00788 rx = (r-REAL(ix-1,dp)*dx)/dx 00789 rx2 = rx*rx 00790 rx3 = rx2*rx 00791 Term = pot0_2(1,ix )*(1.0_dp-3.0_dp*rx2+2.0_dp*rx3) & 00792 +pot0_2(2,ix )*(rx-2.0_dp*rx2+rx3) & 00793 +pot0_2(1,ix+1)*(3.0_dp*rx2-2.0_dp*rx3) & 00794 +pot0_2(2,ix+1)*(-rx2+rx3) 00795 grid2(i,j,k) = grid2(i,j,k) - Term * qt 00796 xs1 = xs1 + dr1 00797 END DO 00798 xs2 = xs2 + dr2 00799 END DO 00800 END DO LoopOnGrid 00801 END DO Atoms 00802 END DO Radius 00803 CALL timestop(handle) 00804 END SUBROUTINE qmmm_elec_with_gaussian_LR 00805 00806 END MODULE qmmm_gpw_energy
1.7.3