|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00011 MODULE qmmm_se_energy 00012 USE atomic_kind_types, ONLY: atomic_kind_type,& 00013 get_atomic_kind 00014 USE cell_types, ONLY: cell_type,& 00015 pbc 00016 USE cp_control_types, ONLY: dft_control_type 00017 USE cp_dbcsr_interface, ONLY: cp_dbcsr_copy,& 00018 cp_dbcsr_get_block_p,& 00019 cp_dbcsr_init,& 00020 cp_dbcsr_set 00021 USE cp_dbcsr_operations, ONLY: cp_dbcsr_allocate_matrix_set 00022 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 00023 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00024 USE cp_output_handling, ONLY: cp_p_file,& 00025 cp_print_key_finished_output,& 00026 cp_print_key_should_output,& 00027 cp_print_key_unit_nr 00028 USE cp_para_types, ONLY: cp_para_env_type 00029 USE f77_blas 00030 USE input_constants, ONLY: & 00031 do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, & 00032 do_method_pdg, do_method_pm3, do_method_pm6, do_method_pnnl, & 00033 do_method_rm1, do_multipole_none, do_qmmm_coulomb, do_qmmm_gauss, & 00034 do_qmmm_none, do_qmmm_swave, use_orb_basis_set 00035 USE kinds, ONLY: dp 00036 USE message_passing, ONLY: mp_sum 00037 USE particle_types, ONLY: particle_type 00038 USE qmmm_types, ONLY: qmmm_env_qm_type,& 00039 qmmm_pot_p_type,& 00040 qmmm_pot_type 00041 USE qmmm_util, ONLY: spherical_cutoff_factor 00042 USE qs_environment_types, ONLY: get_qs_env,& 00043 qs_environment_type,& 00044 set_qs_env 00045 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type 00046 USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists 00047 USE qs_overlap, ONLY: build_overlap_matrix 00048 USE semi_empirical_int_arrays, ONLY: se_orbital_pointer 00049 USE semi_empirical_integrals, ONLY: corecore,& 00050 rotnuc 00051 USE semi_empirical_types, ONLY: get_se_param,& 00052 se_int_control_type,& 00053 se_taper_type,& 00054 semi_empirical_create,& 00055 semi_empirical_release,& 00056 semi_empirical_type,& 00057 setup_se_int_control_type 00058 USE semi_empirical_utils, ONLY: get_se_type,& 00059 se_param_set_default 00060 USE termination, ONLY: stop_program 00061 USE timings, ONLY: timeset,& 00062 timestop 00063 #include "cp_common_uses.h" 00064 00065 IMPLICIT NONE 00066 00067 PRIVATE 00068 00069 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_se_energy' 00070 00071 PUBLIC :: build_se_qmmm_matrix 00072 00073 CONTAINS 00074 00075 ! ***************************************************************************** 00079 SUBROUTINE build_se_qmmm_matrix(qs_env,qmmm_env,mm_particles,mm_cell,para_env,error) 00080 00081 TYPE(qs_environment_type), POINTER :: qs_env 00082 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 00083 TYPE(particle_type), DIMENSION(:), 00084 POINTER :: mm_particles 00085 TYPE(cell_type), POINTER :: mm_cell 00086 TYPE(cp_para_env_type), POINTER :: para_env 00087 TYPE(cp_error_type), INTENT(inout) :: error 00088 00089 CHARACTER(len=*), PARAMETER :: routineN = 'build_se_qmmm_matrix', 00090 routineP = moduleN//':'//routineN 00091 00092 INTEGER :: handle, i, iatom, ikind, 00093 itype, iw, natom, natorb_a, 00094 nkind, stat 00095 INTEGER, DIMENSION(:), POINTER :: list 00096 LOGICAL :: anag, defined, failure, found 00097 REAL(KIND=dp) :: enuclear 00098 REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block_a 00099 TYPE(atomic_kind_type), DIMENSION(:), 00100 POINTER :: atomic_kind_set 00101 TYPE(atomic_kind_type), POINTER :: atomic_kind 00102 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00103 POINTER :: matrix_s 00104 TYPE(cp_logger_type), POINTER :: logger 00105 TYPE(dft_control_type), POINTER :: dft_control 00106 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc 00107 TYPE(se_int_control_type) :: se_int_control 00108 TYPE(se_taper_type), POINTER :: se_taper 00109 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm 00110 00111 failure = .FALSE. 00112 CALL timeset(routineN,handle) 00113 NULLIFY(logger) 00114 logger => cp_error_get_logger(error) 00115 00116 NULLIFY (matrix_s, atomic_kind_set) 00117 NULLIFY (se_kind_a, se_kind_mm, se_taper) 00118 CALL build_qs_neighbor_lists(qs_env,para_env,force_env_section=qs_env%input,error=error) 00119 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, error=error) 00120 CALL build_overlap_matrix(qs_env,matrix_s=matrix_s,& 00121 matrix_name="OVERLAP",& 00122 basis_set_id_a=use_orb_basis_set,& 00123 basis_set_id_b=use_orb_basis_set, & 00124 sab_nl=qs_env%sab_orb,& 00125 error=error) 00126 00127 CALL set_qs_env(qs_env=qs_env,matrix_s=matrix_s,error=error) 00128 CALL get_qs_env(qs_env=qs_env,& 00129 se_taper=se_taper,& 00130 atomic_kind_set=atomic_kind_set,& 00131 ks_qmmm_env=ks_qmmm_env_loc,& 00132 dft_control=dft_control,error=error) 00133 00134 SELECT CASE (dft_control%qs_control%method_id) 00135 CASE (do_method_am1,do_method_rm1,do_method_mndo,do_method_pdg,& 00136 do_method_pm3,do_method_pm6,do_method_mndod, do_method_pnnl) 00137 ! Go on with the calculation.. 00138 CASE DEFAULT 00139 ! Otherwise stop.. 00140 CALL stop_program(routineN,moduleN,__LINE__,& 00141 "Method not available",para_env) 00142 END SELECT 00143 anag = dft_control%qs_control%se_control%analytical_gradients 00144 ! Setup type for SE integral control 00145 CALL setup_se_int_control_type(se_int_control, shortrange=.FALSE.,do_ewald_r3=.FALSE.,& 00146 do_ewald_gks=.FALSE., integral_screening=dft_control%qs_control%se_control%integral_screening,& 00147 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 00148 00149 ! Allocate the core Hamiltonian matrix 00150 CALL cp_dbcsr_allocate_matrix_set(ks_qmmm_env_loc%matrix_h,1,error=error) 00151 ALLOCATE(ks_qmmm_env_loc%matrix_h(1)%matrix, STAT=stat) 00152 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00153 00154 CALL cp_dbcsr_init(ks_qmmm_env_loc%matrix_h(1)%matrix,error) 00155 00156 CALL cp_dbcsr_copy(ks_qmmm_env_loc%matrix_h(1)%matrix,matrix_s(1)%matrix,& 00157 name="QMMM HAMILTONIAN MATRIX",error=error) 00158 CALL cp_dbcsr_set(ks_qmmm_env_loc%matrix_h(1)%matrix,0.0_dp,error) 00159 00160 SELECT CASE(qmmm_env%qmmm_coupl_type) 00161 CASE(do_qmmm_coulomb,do_qmmm_gauss,do_qmmm_swave) 00162 ! Create a fake semi-empirical type to handle the classical atom 00163 CALL semi_empirical_create(se_kind_mm,error) 00164 CALL se_param_set_default(se_kind_mm,0,do_method_pchg,error) 00165 itype = get_se_type(se_kind_mm%typ) 00166 nkind = SIZE(atomic_kind_set) 00167 enuclear = 0.0_dp 00168 Kinds: DO ikind=1,nkind 00169 atomic_kind => atomic_kind_set(ikind) 00170 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00171 natom=natom,& 00172 se_parameter=se_kind_a,& 00173 atom_list=list) 00174 CALL get_se_param(se_kind_a,& 00175 defined=defined,& 00176 natorb=natorb_a) 00177 IF (.NOT.defined .OR. natorb_a < 1) CYCLE 00178 Atoms: DO i = 1, SIZE(list) 00179 iatom = list(i) 00180 ! Give back block 00181 NULLIFY(h_block_a) 00182 CALL cp_dbcsr_get_block_p(matrix=ks_qmmm_env_loc%matrix_h(1)%matrix,& 00183 row=iatom,col=iatom,BLOCK=h_block_a,found=found) 00184 00185 IF (ASSOCIATED(h_block_a)) THEN 00186 h_block_a = 0.0_dp 00187 ! Real QM/MM computation 00188 CALL build_se_qmmm_matrix_low(h_block_a,& 00189 se_kind_a,& 00190 se_kind_mm,& 00191 qmmm_env%Potentials,& 00192 mm_particles,& 00193 qmmm_env%mm_atom_chrg,& 00194 qmmm_env%mm_el_pot_radius,& 00195 qmmm_env%mm_atom_index,& 00196 qmmm_env%num_mm_atoms,& 00197 mm_cell,& 00198 qmmm_env%qm_atom_index(iatom),& 00199 enuclear,& 00200 itype,& 00201 se_taper,& 00202 se_int_control,& 00203 anag,& 00204 qmmm_env%spherical_cutoff,& 00205 error) 00206 ! Possibly added charges 00207 IF (qmmm_env%move_mm_charges.OR.qmmm_env%add_mm_charges) THEN 00208 CALL build_se_qmmm_matrix_low(h_block_a,& 00209 se_kind_a,& 00210 se_kind_mm,& 00211 qmmm_env%added_charges%potentials,& 00212 qmmm_env%added_charges%added_particles,& 00213 qmmm_env%added_charges%mm_atom_chrg,& 00214 qmmm_env%added_charges%mm_el_pot_radius,& 00215 qmmm_env%added_charges%mm_atom_index,& 00216 qmmm_env%added_charges%num_mm_atoms,& 00217 mm_cell,& 00218 qmmm_env%qm_atom_index(iatom),& 00219 enuclear,& 00220 itype,& 00221 se_taper,& 00222 se_int_control,& 00223 anag,& 00224 qmmm_env%spherical_cutoff,& 00225 error) 00226 END IF 00227 END IF 00228 END DO Atoms 00229 END DO Kinds 00230 CALL mp_sum(enuclear,para_env%group) 00231 qs_env%energy%qmmm_nu = enuclear 00232 CALL semi_empirical_release(se_kind_mm,error) 00233 CASE(do_qmmm_none) 00234 ! Zero Matrix 00235 CALL cp_dbcsr_set(ks_qmmm_env_loc%matrix_h(1)%matrix,0.0_dp,error=error) 00236 END SELECT 00237 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00238 qs_env%input,"QMMM%PRINT%QMMM_MATRIX",error=error),cp_p_file)) THEN 00239 iw = cp_print_key_unit_nr(logger,qs_env%input,"QMMM%PRINT%QMMM_MATRIX",& 00240 extension=".Log",error=error) 00241 CALL cp_dbcsr_write_sparse_matrix(ks_qmmm_env_loc%matrix_h(1)%matrix,4,6,qs_env,para_env,& 00242 scale=1.0_dp,output_unit=iw,error=error) 00243 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 00244 "QMMM%PRINT%QMMM_MATRIX", error=error) 00245 END IF 00246 00247 CALL timestop(handle) 00248 00249 END SUBROUTINE build_se_qmmm_matrix 00250 00251 ! ***************************************************************************** 00255 SUBROUTINE build_se_qmmm_matrix_low(h_block_a, se_kind_a, se_kind_mm, potentials,& 00256 mm_particles, mm_charges, mm_el_pot_radius, mm_atom_index, num_mm_atoms,& 00257 mm_cell, IndQM, enuclear, itype, se_taper, se_int_control, anag, & 00258 qmmm_spherical_cutoff, error) 00259 00260 REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block_a 00261 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm 00262 TYPE(qmmm_pot_p_type), DIMENSION(:), 00263 POINTER :: potentials 00264 TYPE(particle_type), DIMENSION(:), 00265 POINTER :: mm_particles 00266 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges, mm_el_pot_radius 00267 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 00268 INTEGER, INTENT(IN) :: num_mm_atoms 00269 TYPE(cell_type), POINTER :: mm_cell 00270 INTEGER, INTENT(IN) :: IndQM 00271 REAL(KIND=dp), INTENT(INOUT) :: enuclear 00272 INTEGER, INTENT(IN) :: itype 00273 TYPE(se_taper_type), POINTER :: se_taper 00274 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00275 LOGICAL, INTENT(IN) :: anag 00276 REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2) 00277 TYPE(cp_error_type), INTENT(inout) :: error 00278 00279 CHARACTER(len=*), PARAMETER :: routineN = 'build_se_qmmm_matrix_low', 00280 routineP = moduleN//':'//routineN 00281 00282 INTEGER :: handle, i1, i1L, i2, Imm, 00283 Imp, IndMM, Ipot, j1, j1L 00284 REAL(KIND=dp) :: enuc, rt1, rt2, rt3, 00285 sph_chrg_factor 00286 REAL(KIND=dp), DIMENSION(3) :: r_pbc, rij 00287 REAL(KIND=dp), DIMENSION(45) :: e1b 00288 TYPE(qmmm_pot_type), POINTER :: Pot 00289 00290 CALL timeset(routineN,handle) 00291 ! Loop Over MM atoms 00292 ! Loop over Pot stores atoms with the same charge 00293 MainLoopPot: DO Ipot = 1, SIZE(Potentials) 00294 Pot => Potentials(Ipot)%Pot 00295 ! Loop over atoms belonging to this type 00296 LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index) 00297 Imm = Pot%mm_atom_index(Imp) 00298 IndMM = mm_atom_index(Imm) 00299 r_pbc=pbc(mm_particles(IndMM)%r-mm_particles(IndQM)%r, mm_cell) 00300 rt1= r_pbc(1) 00301 rt2= r_pbc(2) 00302 rt3= r_pbc(3) 00303 rij = (/rt1,rt2,rt3/) 00304 se_kind_mm%zeff = mm_charges(Imm) 00305 ! Computes the screening factor for the spherical cutoff (if defined) 00306 IF (qmmm_spherical_cutoff(1)>0.0_dp) THEN 00307 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor, error) 00308 se_kind_mm%zeff = se_kind_mm%zeff * sph_chrg_factor 00309 END IF 00310 IF (ABS(se_kind_mm%zeff)<=EPSILON(0.0_dp)) CYCLE 00311 CALL rotnuc (se_kind_a, se_kind_mm, rij, itype=itype, e1b=e1b, anag=anag,& 00312 se_int_control=se_int_control, se_taper=se_taper, error=error) 00313 CALL corecore(se_kind_a, se_kind_mm, rij, itype=itype, enuc=enuc, anag=anag,& 00314 se_int_control=se_int_control, se_taper=se_taper, error=error) 00315 enuclear = enuclear + enuc 00316 ! Contribution to the iatom block 00317 ! Computation of the QMMM core matrix 00318 i2 = 0 00319 DO i1L = 1, se_kind_a%natorb 00320 i1 = se_orbital_pointer(i1L) 00321 DO j1L = 1, i1L-1 00322 j1 = se_orbital_pointer(j1L) 00323 i2 = i2 + 1 00324 h_block_a(i1,j1) = h_block_a(i1,j1) + e1b(i2) 00325 h_block_a(j1,i1) = h_block_a(i1,j1) 00326 END DO 00327 j1 = se_orbital_pointer(j1L) 00328 i2 = i2 + 1 00329 h_block_a(i1,j1) = h_block_a(i1,j1) + e1b(i2) 00330 END DO 00331 END DO LoopMM 00332 END DO MainLoopPot 00333 CALL timestop(handle) 00334 END SUBROUTINE build_se_qmmm_matrix_low 00335 00336 END MODULE qmmm_se_energy
1.7.3