|
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 ! ***************************************************************************** 00017 MODULE se_fock_matrix_coulomb_mpi 00018 USE atomic_kind_types, ONLY: atomic_kind_type,& 00019 get_atomic_kind,& 00020 get_atomic_kind_set 00021 USE cell_types, ONLY: cell_type 00022 USE cp_control_types, ONLY: dft_control_type,& 00023 semi_empirical_control_type 00024 USE cp_dbcsr_interface, ONLY: cp_dbcsr_add,& 00025 cp_dbcsr_distribute,& 00026 cp_dbcsr_get_block_diag,& 00027 cp_dbcsr_get_block_p,& 00028 cp_dbcsr_init,& 00029 cp_dbcsr_replicate_all,& 00030 cp_dbcsr_set,& 00031 cp_dbcsr_sum_replicated 00032 USE cp_dbcsr_operations, ONLY: cp_dbcsr_allocate_matrix_set,& 00033 cp_dbcsr_deallocate_matrix_set 00034 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00035 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00036 cp_print_key_unit_nr 00037 USE cp_para_types, ONLY: cp_para_env_type 00038 USE distribution_1d_types, ONLY: distribution_1d_type 00039 USE ewald_environment_types, ONLY: ewald_env_get,& 00040 ewald_environment_type 00041 USE ewald_pw_types, ONLY: ewald_pw_get,& 00042 ewald_pw_type 00043 USE ewalds_multipole, ONLY: ewald_multipole_evaluate 00044 USE f77_blas 00045 USE fist_neighbor_list_control, ONLY: list_control 00046 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type 00047 USE input_constants, ONLY: & 00048 do_ewald_ewald, do_method_am1, do_method_mndo, do_method_mndod, & 00049 do_method_pdg, do_method_pm3, do_method_pm6, do_method_pnnl, & 00050 do_method_rm1, do_multipole_charge, do_multipole_dipole, & 00051 do_multipole_none, do_multipole_quadrupole, do_se_IS_slater 00052 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00053 section_vals_type 00054 USE kinds, ONLY: dp 00055 USE message_passing, ONLY: mp_sum 00056 USE particle_types, ONLY: particle_type 00057 USE qs_energy_types, ONLY: qs_energy_type 00058 USE qs_environment_types, ONLY: get_qs_env,& 00059 qs_environment_type 00060 USE qs_force_types, ONLY: qs_force_type 00061 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 00062 neighbor_list_iterate,& 00063 neighbor_list_iterator_create,& 00064 neighbor_list_iterator_p_type,& 00065 neighbor_list_iterator_release,& 00066 neighbor_list_set_p_type 00067 USE se_fock_matrix_integrals, ONLY: & 00068 dfock2C, dfock2C_r3, dfock2_1el, dfock2_1el_r3, fock2C, fock2C_ew, & 00069 fock2C_r3, fock2_1el, fock2_1el_ew, fock2_1el_r3, & 00070 se_coulomb_ij_interaction 00071 USE semi_empirical_int_arrays, ONLY: rij_threshold,& 00072 se_orbital_pointer 00073 USE semi_empirical_integrals, ONLY: corecore_el,& 00074 dcorecore_el 00075 USE semi_empirical_mpole_methods, ONLY: quadrupole_sph_to_cart 00076 USE semi_empirical_mpole_types, ONLY: nddo_mpole_type,& 00077 semi_empirical_mpole_type 00078 USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type 00079 USE semi_empirical_types, ONLY: get_se_param,& 00080 se_int_control_type,& 00081 se_taper_type,& 00082 semi_empirical_p_type,& 00083 semi_empirical_type,& 00084 setup_se_int_control_type 00085 USE semi_empirical_utils, ONLY: finalize_se_taper,& 00086 get_se_type,& 00087 initialize_se_taper 00088 USE timings, ONLY: timeset,& 00089 timestop 00090 USE virial_methods, ONLY: virial_pair_force 00091 USE virial_types, ONLY: virial_type 00092 #include "cp_common_uses.h" 00093 00094 IMPLICIT NONE 00095 PRIVATE 00096 00097 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_coulomb_mpi' 00098 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 00099 00100 PUBLIC :: build_fock_matrix_coulomb, build_fock_matrix_coulomb_lr,& 00101 build_fock_matrix_coul_lr_r3 00102 00103 CONTAINS 00104 00105 ! ***************************************************************************** 00109 SUBROUTINE build_fock_matrix_coulomb (qs_env,ks_matrix,matrix_p,energy,calculate_forces,& 00110 store_int_env,error) 00111 00112 TYPE(qs_environment_type), POINTER :: qs_env 00113 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00114 POINTER :: ks_matrix, matrix_p 00115 TYPE(qs_energy_type), POINTER :: energy 00116 LOGICAL, INTENT(in) :: calculate_forces 00117 TYPE(semi_empirical_si_type), POINTER :: store_int_env 00118 TYPE(cp_error_type), INTENT(inout) :: error 00119 00120 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coulomb', 00121 routineP = moduleN//':'//routineN 00122 00123 INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, ispin, itype, 00124 jatom, jkind, natom, natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, 00125 nspins, stat 00126 INTEGER, ALLOCATABLE, DIMENSION(:) :: se_natorb 00127 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 00128 LOGICAL :: anag, atener, check, defined, 00129 failure, found, switch, 00130 use_virial 00131 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 00132 REAL(KIND=dp) :: delta, dr1, ecore2, ecores 00133 REAL(KIND=dp), DIMENSION(2) :: ecab 00134 REAL(KIND=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b 00135 REAL(KIND=dp), DIMENSION(3) :: force_ab, rij 00136 REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b 00137 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, 00138 ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, pb_block_a, pb_block_b 00139 TYPE(atomic_kind_type), DIMENSION(:), 00140 POINTER :: atomic_kind_set 00141 TYPE(atomic_kind_type), POINTER :: atomic_kind 00142 TYPE(cell_type), POINTER :: cell 00143 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00144 POINTER :: diagmat_ks, diagmat_p 00145 TYPE(cp_para_env_type), POINTER :: para_env 00146 TYPE(dft_control_type), POINTER :: dft_control 00147 TYPE(ewald_environment_type), POINTER :: ewald_env 00148 TYPE(ewald_pw_type), POINTER :: ewald_pw 00149 TYPE(neighbor_list_iterator_p_type), 00150 DIMENSION(:), POINTER :: nl_iterator 00151 TYPE(neighbor_list_set_p_type), 00152 DIMENSION(:), POINTER :: sab_se 00153 TYPE(particle_type), DIMENSION(:), 00154 POINTER :: particle_set 00155 TYPE(qs_force_type), DIMENSION(:), 00156 POINTER :: force 00157 TYPE(se_int_control_type) :: se_int_control 00158 TYPE(se_taper_type), POINTER :: se_taper 00159 TYPE(semi_empirical_control_type), 00160 POINTER :: se_control 00161 TYPE(semi_empirical_p_type), 00162 DIMENSION(:), POINTER :: se_kind_list 00163 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 00164 TYPE(virial_type), POINTER :: virial 00165 00166 failure=.FALSE. 00167 CALL timeset(routineN,handle) 00168 00169 IF ( .NOT. failure ) THEN 00170 NULLIFY(dft_control,cell,force,particle_set,diagmat_ks,diagmat_p,& 00171 se_control,se_taper) 00172 00173 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper,& 00174 para_env=para_env, sab_se=sab_se, atomic_kind_set=atomic_kind_set, & 00175 particle_set=particle_set, virial=virial, error=error) 00176 00177 ! Parameters 00178 CALL initialize_se_taper(se_taper,coulomb=.TRUE.,error=error) 00179 se_control => dft_control%qs_control%se_control 00180 anag = se_control%analytical_gradients 00181 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) 00182 atener = qs_env%atprop%energy 00183 00184 CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3,& 00185 do_ewald_gks=se_control%do_ewald_gks, integral_screening=se_control%integral_screening,& 00186 shortrange=(se_control%do_ewald.OR.se_control%do_ewald_gks),& 00187 max_multipole=se_control%max_multipole,pc_coulomb_int=.FALSE.) 00188 IF(se_control%do_ewald_gks) THEN 00189 CALL get_qs_env(qs_env=qs_env,ewald_env=ewald_env,ewald_pw=ewald_pw, error=error) 00190 CALL ewald_env_get (ewald_env, alpha=se_int_control%ewald_gks%alpha, error=error) 00191 CALL ewald_pw_get (ewald_pw, pw_big_pool=se_int_control%ewald_gks%pw_pool, & 00192 dg=se_int_control%ewald_gks%dg) 00193 END IF 00194 00195 nspins=dft_control%nspins 00196 CPPrecondition(ASSOCIATED(matrix_p),cp_failure_level,routineP,error,failure) 00197 CPPrecondition(SIZE(ks_matrix)>0,cp_failure_level,routineP,error,failure) 00198 00199 nkind = SIZE(atomic_kind_set) 00200 IF(calculate_forces) THEN 00201 CALL get_qs_env(qs_env=qs_env, force=force, error=error) 00202 natom = SIZE (particle_set) 00203 delta = se_control%delta 00204 ! Allocate atom index for kind 00205 ALLOCATE (atom_of_kind(natom),STAT=stat) 00206 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00207 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,atom_of_kind=atom_of_kind) 00208 END IF 00209 00210 CALL cp_dbcsr_allocate_matrix_set(diagmat_ks,nspins,error=error) 00211 CALL cp_dbcsr_allocate_matrix_set(diagmat_p,nspins,error=error) 00212 00213 DO ispin=1,nspins 00214 ! Allocate diagonal block matrices 00215 ALLOCATE(diagmat_p(ispin)%matrix,diagmat_ks(ispin)%matrix)!sm->dbcsr 00216 CALL cp_dbcsr_init(diagmat_p(ispin)%matrix, error=error) 00217 CALL cp_dbcsr_init(diagmat_ks(ispin)%matrix, error=error) 00218 CALL cp_dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix, error=error) 00219 CALL cp_dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, error=error) 00220 CALL cp_dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp, error=error) 00221 CALL cp_dbcsr_replicate_all(diagmat_p(ispin)%matrix, error=error) 00222 CALL cp_dbcsr_replicate_all(diagmat_ks(ispin)%matrix, error=error) 00223 END DO 00224 00225 ecore2 = 0.0_dp 00226 itype = get_se_type(dft_control%qs_control%method_id) 00227 ALLOCATE (se_defined(nkind),se_kind_list(nkind),se_natorb(nkind),STAT=stat) 00228 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00229 DO ikind=1,nkind 00230 atomic_kind => atomic_kind_set(ikind) 00231 CALL get_atomic_kind(atomic_kind=atomic_kind,se_parameter=se_kind_a) 00232 se_kind_list(ikind)%se_param => se_kind_a 00233 CALL get_se_param(se_kind_a,defined=defined,natorb=natorb_a) 00234 se_defined(ikind) = (defined .AND. natorb_a >= 1) 00235 se_natorb(ikind) = natorb_a 00236 END DO 00237 00238 CALL neighbor_list_iterator_create(nl_iterator,sab_se) 00239 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 00240 CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,iatom=iatom,jatom=jatom,inode=inode,r=rij) 00241 IF (.NOT.se_defined(ikind)) CYCLE 00242 IF (.NOT.se_defined(jkind)) CYCLE 00243 se_kind_a => se_kind_list(ikind)%se_param 00244 se_kind_b => se_kind_list(jkind)%se_param 00245 natorb_a = se_natorb(ikind) 00246 natorb_b = se_natorb(jkind) 00247 natorb_a2 = natorb_a**2 00248 natorb_b2 = natorb_b**2 00249 00250 IF (inode==1) THEN 00251 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(1)%matrix,& 00252 row=iatom,col=iatom,BLOCK=pa_block_a,found=found) 00253 CPPostcondition(ASSOCIATED(pa_block_a),cp_failure_level,routineP,error,failure) 00254 check = (SIZE(pa_block_a,1)==natorb_a).AND.(SIZE(pa_block_a,2)==natorb_a) 00255 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00256 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix,& 00257 row=iatom,col=iatom,BLOCK=ksa_block_a,found=found) 00258 CPPostcondition(ASSOCIATED(ksa_block_a),cp_failure_level,routineP,error,failure) 00259 p_block_tot_a(1:natorb_a,1:natorb_a) = 2.0_dp * pa_block_a 00260 pa_a(1:natorb_a2) = RESHAPE(pa_block_a,(/natorb_a2/)) 00261 IF ( nspins == 2 ) THEN 00262 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(2)%matrix,& 00263 row=iatom,col=iatom,BLOCK=pa_block_b,found=found) 00264 CPPostcondition(ASSOCIATED(pa_block_b),cp_failure_level,routineP,error,failure) 00265 check = (SIZE(pa_block_b,1)==natorb_a).AND.(SIZE(pa_block_b,2)==natorb_a) 00266 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00267 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix,& 00268 row=iatom,col=iatom,BLOCK=ksa_block_b,found=found) 00269 CPPostcondition(ASSOCIATED(ksa_block_b),cp_failure_level,routineP,error,failure) 00270 p_block_tot_a(1:natorb_a,1:natorb_a) = pa_block_a + pa_block_b 00271 pa_b(1:natorb_a2) = RESHAPE(pa_block_b,(/natorb_a2/)) 00272 END IF 00273 00274 ENDIF 00275 00276 dr1 = DOT_PRODUCT(rij,rij) 00277 IF ( dr1 > rij_threshold ) THEN 00278 ! Determine the order of the atoms, and in case switch them.. 00279 IF (iatom <= jatom) THEN 00280 switch = .FALSE. 00281 ELSE 00282 switch = .TRUE. 00283 END IF 00284 ! Retrieve blocks for KS and P 00285 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(1)%matrix,& 00286 row=jatom,col=jatom,BLOCK=pb_block_a,found=found) 00287 CPPostcondition(ASSOCIATED(pb_block_a),cp_failure_level,routineP,error,failure) 00288 check = (SIZE(pb_block_a,1)==natorb_b).AND.(SIZE(pb_block_a,2)==natorb_b) 00289 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00290 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix,& 00291 row=jatom,col=jatom,BLOCK=ksb_block_a,found=found) 00292 CPPostcondition(ASSOCIATED(ksb_block_a),cp_failure_level,routineP,error,failure) 00293 p_block_tot_b(1:natorb_b,1:natorb_b) = 2.0_dp * pb_block_a 00294 pb_a(1:natorb_b2) = RESHAPE(pb_block_a,(/natorb_b2/)) 00295 ! Handle more than one configuration 00296 IF ( nspins == 2 ) THEN 00297 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(2)%matrix,& 00298 row=jatom,col=jatom,BLOCK=pb_block_b,found=found) 00299 CPPostcondition(ASSOCIATED(pb_block_b),cp_failure_level,routineP,error,failure) 00300 check = (SIZE(pb_block_b,1)==natorb_b).AND.(SIZE(pb_block_b,2)==natorb_b) 00301 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00302 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix,& 00303 row=jatom,col=jatom,BLOCK=ksb_block_b,found=found) 00304 CPPostcondition(ASSOCIATED(ksb_block_b),cp_failure_level,routineP,error,failure) 00305 check = (SIZE(pb_block_a,1)==SIZE(pb_block_b,1)).AND.(SIZE(pb_block_a,2)==SIZE(pb_block_b,2)) 00306 CPPostcondition(check,cp_failure_level,routineP,error,failure) 00307 p_block_tot_b(1:natorb_b,1:natorb_b) = pb_block_a + pb_block_b 00308 pb_b(1:natorb_b2) = RESHAPE(pb_block_b,(/natorb_b2/)) 00309 END IF 00310 00311 SELECT CASE (dft_control%qs_control%method_id) 00312 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pdg,& 00313 do_method_rm1, do_method_mndod, do_method_pnnl) 00314 00315 ! Two-centers One-electron terms 00316 IF ( nspins == 1 ) THEN 00317 ecab = 0._dp 00318 CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_a, ksb_block_a,& 00319 pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control,& 00320 se_taper=se_taper, store_int_env=store_int_env, error=error) 00321 ecore2 = ecore2 + ecab(1) + ecab(2) 00322 ELSE IF ( nspins == 2 ) THEN 00323 ecab = 0._dp 00324 CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_a, ksb_block_a,& 00325 pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag,& 00326 se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env,& 00327 error=error) 00328 CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_b, ksb_block_b,& 00329 pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control,& 00330 se_taper=se_taper, store_int_env=store_int_env, error=error) 00331 ecore2 = ecore2 + ecab(1) + ecab(2) 00332 END IF 00333 IF (atener) THEN 00334 qs_env%atprop%atecoul(iatom) = qs_env%atprop%atecoul(iatom) + ecab(1) 00335 qs_env%atprop%atecoul(jatom) = qs_env%atprop%atecoul(jatom) + ecab(2) 00336 END IF 00337 ! Coulomb Terms 00338 IF ( nspins == 1 ) THEN 00339 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,& 00340 ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 00341 store_int_env=store_int_env, error=error) 00342 ELSE IF ( nspins == 2 ) THEN 00343 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,& 00344 ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 00345 store_int_env=store_int_env, error=error) 00346 00347 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b,& 00348 ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 00349 store_int_env=store_int_env, error=error) 00350 END IF 00351 00352 IF(calculate_forces) THEN 00353 atom_a = atom_of_kind(iatom) 00354 atom_b = atom_of_kind(jatom) 00355 00356 ! Derivatives of the Two-centre One-electron terms 00357 force_ab = 0.0_dp 00358 IF ( nspins == 1 ) THEN 00359 CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_a, pb_a, itype=itype, anag=anag,& 00360 se_int_control=se_int_control, se_taper=se_taper, force=force_ab,& 00361 delta=delta, error=error) 00362 ELSE IF ( nspins == 2 ) THEN 00363 CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_block_a, pb_block_a, itype=itype, anag=anag,& 00364 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 00365 error=error) 00366 CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_b, pb_b, itype=itype, anag=anag,& 00367 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 00368 error=error) 00369 END IF 00370 IF (use_virial) THEN 00371 CALL virial_pair_force (virial%pv_virial, -1.0_dp, force_ab, rij, error) 00372 END IF 00373 00374 ! Sum up force components 00375 force(ikind)%all_potential(1,atom_a) = force(ikind)%all_potential(1,atom_a) - force_ab(1) 00376 force(jkind)%all_potential(1,atom_b) = force(jkind)%all_potential(1,atom_b) + force_ab(1) 00377 00378 force(ikind)%all_potential(2,atom_a) = force(ikind)%all_potential(2,atom_a) - force_ab(2) 00379 force(jkind)%all_potential(2,atom_b) = force(jkind)%all_potential(2,atom_b) + force_ab(2) 00380 00381 force(ikind)%all_potential(3,atom_a) = force(ikind)%all_potential(3,atom_a) - force_ab(3) 00382 force(jkind)%all_potential(3,atom_b) = force(jkind)%all_potential(3,atom_b) + force_ab(3) 00383 00384 ! Derivatives of the Coulomb Terms 00385 force_ab = 0._dp 00386 IF ( nspins == 1 ) THEN 00387 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp,& 00388 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 00389 error=error) 00390 ELSE IF ( nspins == 2 ) THEN 00391 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,& 00392 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 00393 error=error) 00394 00395 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,& 00396 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 00397 error=error) 00398 END IF 00399 IF ( switch ) THEN 00400 force_ab(1) = -force_ab(1) 00401 force_ab(2) = -force_ab(2) 00402 force_ab(3) = -force_ab(3) 00403 END IF 00404 IF (use_virial) THEN 00405 CALL virial_pair_force ( virial%pv_virial, -1.0_dp, force_ab, rij, error) 00406 END IF 00407 ! Sum up force components 00408 force(ikind)%rho_elec(1,atom_a) = force(ikind)%rho_elec(1,atom_a) - force_ab(1) 00409 force(jkind)%rho_elec(1,atom_b) = force(jkind)%rho_elec(1,atom_b) + force_ab(1) 00410 00411 force(ikind)%rho_elec(2,atom_a) = force(ikind)%rho_elec(2,atom_a) - force_ab(2) 00412 force(jkind)%rho_elec(2,atom_b) = force(jkind)%rho_elec(2,atom_b) + force_ab(2) 00413 00414 force(ikind)%rho_elec(3,atom_a) = force(ikind)%rho_elec(3,atom_a) - force_ab(3) 00415 force(jkind)%rho_elec(3,atom_b) = force(jkind)%rho_elec(3,atom_b) + force_ab(3) 00416 END IF 00417 CASE DEFAULT 00418 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00419 END SELECT 00420 ELSE 00421 IF ( se_int_control%do_ewald_gks ) THEN 00422 CPPostcondition(iatom==jatom,cp_failure_level,routineP,error,failure) 00423 ! Two-centers One-electron terms 00424 ecores=0._dp 00425 IF ( nspins == 1 ) THEN 00426 CALL fock2_1el_ew (se_kind_a,rij, ksa_block_a, pa_a, & 00427 ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control,& 00428 se_taper=se_taper, store_int_env=store_int_env, error=error) 00429 ELSE IF ( nspins == 2 ) THEN 00430 CALL fock2_1el_ew (se_kind_a,rij, ksa_block_a, pa_block_a,& 00431 ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, & 00432 se_taper=se_taper, store_int_env=store_int_env, error=error) 00433 CALL fock2_1el_ew (se_kind_a,rij, ksa_block_b, pa_b,& 00434 ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control,& 00435 se_taper=se_taper, store_int_env=store_int_env, error=error) 00436 END IF 00437 ecore2=ecore2+ecores 00438 IF (atener) THEN 00439 qs_env%atprop%atecoul(iatom) = qs_env%atprop%atecoul(iatom) + ecores 00440 END IF 00441 ! Coulomb Terms 00442 IF ( nspins == 1 ) THEN 00443 CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a,& 00444 factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 00445 store_int_env=store_int_env, error=error) 00446 ELSE IF ( nspins == 2 ) THEN 00447 CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a,& 00448 factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 00449 store_int_env=store_int_env, error=error) 00450 CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_b,& 00451 factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 00452 store_int_env=store_int_env, error=error) 00453 END IF 00454 END IF 00455 END IF 00456 END DO 00457 CALL neighbor_list_iterator_release(nl_iterator) 00458 00459 DEALLOCATE(se_kind_list,se_defined,se_natorb,stat=stat) 00460 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00461 00462 DO ispin = 1, nspins 00463 CALL cp_dbcsr_sum_replicated( diagmat_ks(ispin)%matrix, error=error) 00464 CALL cp_dbcsr_distribute(diagmat_ks(ispin)%matrix, error=error) 00465 CALL cp_dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix,& 00466 1.0_dp,1.0_dp,error=error) 00467 END DO 00468 CALL cp_dbcsr_deallocate_matrix_set ( diagmat_p, error=error ) 00469 CALL cp_dbcsr_deallocate_matrix_set ( diagmat_ks, error=error ) 00470 00471 IF (calculate_forces) THEN 00472 DEALLOCATE(atom_of_kind,stat=stat) 00473 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00474 END IF 00475 00476 ! Two-centers one-electron terms 00477 CALL mp_sum(ecore2,para_env%group) 00478 energy%hartree = ecore2 - energy%core 00479 ! WRITE ( *, * ) 'IN SE_F_COUL', ecore2, energy%core 00480 00481 END IF 00482 CALL finalize_se_taper(se_taper,error=error) 00483 00484 CALL timestop(handle) 00485 END SUBROUTINE build_fock_matrix_coulomb 00486 00487 ! ***************************************************************************** 00492 SUBROUTINE build_fock_matrix_coulomb_lr (qs_env, ks_matrix, matrix_p, energy,& 00493 calculate_forces, store_int_env, error) 00494 00495 TYPE(qs_environment_type), POINTER :: qs_env 00496 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00497 POINTER :: ks_matrix, matrix_p 00498 TYPE(qs_energy_type), POINTER :: energy 00499 LOGICAL, INTENT(in) :: calculate_forces 00500 TYPE(semi_empirical_si_type), POINTER :: store_int_env 00501 TYPE(cp_error_type), INTENT(inout) :: error 00502 00503 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coulomb_lr', 00504 routineP = moduleN//':'//routineN 00505 00506 INTEGER :: atom_a, ewald_type, forces_g_size, handle, iatom, ikind, 00507 ilist, indi, indj, ispin, itype, iw, jint, natoms, natorb_a, nkind, 00508 nlocal_particles, node, nparticle_local, nspins, size_1c_int, stat 00509 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 00510 LOGICAL :: anag, atener, defined, 00511 failure, found, use_virial 00512 LOGICAL, DIMENSION(3) :: task 00513 REAL(KIND=dp) :: e_neut, e_self, energy_glob, 00514 energy_local, enuc, fac, tmp 00515 REAL(KIND=dp), ALLOCATABLE, 00516 DIMENSION(:, :) :: forces_g, forces_r 00517 REAL(KIND=dp), DIMENSION(3) :: force_a 00518 REAL(KIND=dp), DIMENSION(3, 3) :: pv_glob, pv_local, qcart 00519 REAL(KIND=dp), DIMENSION(5) :: qsph 00520 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, pa_block_a 00521 TYPE(atomic_kind_type), DIMENSION(:), 00522 POINTER :: atomic_kind_set 00523 TYPE(atomic_kind_type), POINTER :: atomic_kind 00524 TYPE(cell_type), POINTER :: cell 00525 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00526 POINTER :: diagmat_ks, diagmat_p 00527 TYPE(cp_logger_type), POINTER :: logger 00528 TYPE(cp_para_env_type), POINTER :: para_env 00529 TYPE(dft_control_type), POINTER :: dft_control 00530 TYPE(distribution_1d_type), POINTER :: local_particles 00531 TYPE(ewald_environment_type), POINTER :: ewald_env 00532 TYPE(ewald_pw_type), POINTER :: ewald_pw 00533 TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env 00534 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 00535 TYPE(particle_type), DIMENSION(:), 00536 POINTER :: particle_set 00537 TYPE(qs_force_type), DIMENSION(:), 00538 POINTER :: force 00539 TYPE(section_vals_type), POINTER :: se_section 00540 TYPE(semi_empirical_control_type), 00541 POINTER :: se_control 00542 TYPE(semi_empirical_mpole_type), POINTER :: mpole 00543 TYPE(semi_empirical_type), POINTER :: se_kind_a 00544 TYPE(virial_type), POINTER :: virial 00545 00546 failure=.FALSE. 00547 CALL timeset(routineN,handle) 00548 00549 IF ( .NOT. failure ) THEN 00550 NULLIFY(dft_control, cell, force,particle_set, diagmat_ks, diagmat_p,local_particles,& 00551 se_control, ewald_env, ewald_pw, se_nddo_mpole, se_nonbond_env, se_section, mpole,& 00552 logger) 00553 00554 logger => cp_error_get_logger(error) 00555 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env,& 00556 atomic_kind_set=atomic_kind_set, particle_set=particle_set, ewald_env=ewald_env,& 00557 local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole,& 00558 se_nonbond_env=se_nonbond_env, virial=virial, & 00559 error=error) 00560 00561 nlocal_particles = SUM(local_particles%n_el(:)) 00562 natoms = SIZE(particle_set) 00563 CALL ewald_env_get (ewald_env, ewald_type=ewald_type, error=error) 00564 SELECT CASE(ewald_type) 00565 CASE (do_ewald_ewald) 00566 forces_g_size= nlocal_particles 00567 CASE DEFAULT 00568 CALL cp_unimplemented_error(fromWhere=routineP, & 00569 message="Periodic SE implemented only for standard EWALD sums.", & 00570 error=error, error_level=cp_failure_level) 00571 END SELECT 00572 00573 ! Parameters 00574 se_section => section_vals_get_subs_vals(qs_env%input,"DFT%QS%SE",error=error) 00575 se_control => dft_control%qs_control%se_control 00576 anag = se_control%analytical_gradients 00577 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer).AND.calculate_forces 00578 atener = qs_env%atprop%energy 00579 00580 nspins=dft_control%nspins 00581 CPPrecondition(ASSOCIATED(matrix_p),cp_failure_level,routineP,error,failure) 00582 CPPrecondition(SIZE(ks_matrix)>0,cp_failure_level,routineP,error,failure) 00583 00584 CALL cp_dbcsr_allocate_matrix_set(diagmat_ks,nspins,error=error) 00585 CALL cp_dbcsr_allocate_matrix_set(diagmat_p,nspins,error=error) 00586 00587 nkind = SIZE(atomic_kind_set) 00588 DO ispin=1,nspins 00589 ! Allocate diagonal block matrices 00590 ALLOCATE(diagmat_p(ispin)%matrix,diagmat_ks(ispin)%matrix)!sm->dbcsr 00591 CALL cp_dbcsr_init(diagmat_p(ispin)%matrix, error=error) 00592 CALL cp_dbcsr_init(diagmat_ks(ispin)%matrix, error=error) 00593 CALL cp_dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix, error=error) 00594 CALL cp_dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, error=error) 00595 CALL cp_dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp, error=error) 00596 CALL cp_dbcsr_replicate_all(diagmat_p(ispin)%matrix, error=error) 00597 CALL cp_dbcsr_replicate_all(diagmat_ks(ispin)%matrix, error=error) 00598 END DO 00599 00600 ! Check for implemented SE methods 00601 SELECT CASE (dft_control%qs_control%method_id) 00602 CASE(do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pdg,& 00603 do_method_rm1, do_method_mndod, do_method_pnnl) 00604 itype = get_se_type(dft_control%qs_control%method_id) 00605 CASE DEFAULT 00606 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00607 END SELECT 00608 00609 ! Zero arrays and possibly build neighbor lists 00610 energy_local = 0.0_dp 00611 energy_glob = 0.0_dp 00612 e_neut = 0.0_dp 00613 e_self = 0.0_dp 00614 task = .FALSE. 00615 SELECT CASE(se_control%max_multipole) 00616 CASE (do_multipole_none) 00617 ! Do Nothing 00618 CASE(do_multipole_charge) 00619 task(1) = .TRUE. 00620 CASE(do_multipole_dipole) 00621 task = .TRUE. 00622 task(3) = .FALSE. 00623 CASE(do_multipole_quadrupole) 00624 task = .TRUE. 00625 CASE DEFAULT 00626 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00627 END SELECT 00628 00629 ! Build-up neighbor lists for real-space part of Ewald multipoles 00630 CALL list_control ( atomic_kind_set, particle_set, local_particles, & 00631 cell, se_nonbond_env, para_env, se_section, error=error) 00632 00633 enuc = 0.0_dp 00634 energy%core_overlap = 0.0_dp 00635 se_nddo_mpole%charge = 0.0_dp 00636 se_nddo_mpole%dipole = 0.0_dp 00637 se_nddo_mpole%quadrupole = 0.0_dp 00638 00639 DO ispin = 1, nspins 00640 ! Compute the NDDO mpole expansion 00641 DO ikind=1,nkind 00642 atomic_kind => atomic_kind_set(ikind) 00643 CALL get_atomic_kind(atomic_kind=atomic_kind,se_parameter=se_kind_a) 00644 CALL get_se_param(se_kind_a,defined=defined,natorb=natorb_a) 00645 00646 IF (.NOT.defined .OR. natorb_a < 1) CYCLE 00647 00648 nparticle_local = local_particles%n_el(ikind) 00649 DO ilist=1, nparticle_local 00650 iatom = local_particles%list(ikind)%array(ilist) 00651 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(ispin)%matrix,& 00652 row=iatom,col=iatom,BLOCK=pa_block_a,found=found) 00653 CPPostcondition(ASSOCIATED(pa_block_a),cp_failure_level,routineP,error,failure) 00654 ! Nuclei 00655 IF (task(1).AND.ispin==1) se_nddo_mpole%charge(iatom) = se_kind_a%zeff 00656 ! Electrons 00657 size_1c_int = SIZE(se_kind_a%w_mpole) 00658 DO jint = 1, size_1c_int 00659 mpole => se_kind_a%w_mpole(jint)%mpole 00660 indi = se_orbital_pointer(mpole%indi) 00661 indj = se_orbital_pointer(mpole%indj) 00662 fac = 1.0_dp 00663 IF (indi/=indj) fac = 2.0_dp 00664 00665 ! Charge 00666 IF (mpole%task(1).AND.task(1)) THEN 00667 se_nddo_mpole%charge(iatom) = se_nddo_mpole%charge(iatom) +& 00668 fac*pa_block_a(indi,indj)*mpole%c 00669 END IF 00670 00671 ! Dipole 00672 IF (mpole%task(2).AND.task(2)) THEN 00673 se_nddo_mpole%dipole(:,iatom) = se_nddo_mpole%dipole(:,iatom) +& 00674 fac*pa_block_a(indi,indj)*mpole%d(:) 00675 END IF 00676 00677 ! Quadrupole 00678 IF (mpole%task(3).AND.task(3)) THEN 00679 qsph = fac*mpole%qs * pa_block_a(indi,indj) 00680 CALL quadrupole_sph_to_cart(qcart, qsph, error) 00681 se_nddo_mpole%quadrupole(:,:,iatom) = se_nddo_mpole%quadrupole(:,:,iatom) +& 00682 qcart 00683 END IF 00684 END DO 00685 ! Print some info about charge, dipole and quadrupole (debug purpose only) 00686 IF (debug_this_module) THEN 00687 WRITE(*,'(I5,F12.6,5X,3F12.6,5X,9F12.6)')iatom, se_nddo_mpole%charge(iatom),& 00688 se_nddo_mpole%dipole(:,iatom), se_nddo_mpole%quadrupole(:,:,iatom) 00689 END IF 00690 END DO 00691 END DO 00692 END DO 00693 CALL mp_sum(se_nddo_mpole%charge, para_env%group) 00694 CALL mp_sum(se_nddo_mpole%dipole, para_env%group) 00695 CALL mp_sum(se_nddo_mpole%quadrupole, para_env%group) 00696 00697 ! Initialize for virial 00698 IF (use_virial) THEN 00699 pv_glob = 0.0_dp 00700 pv_local = 0.0_dp 00701 END IF 00702 00703 ! Ewald Multipoles Sum 00704 iw = cp_print_key_unit_nr(logger,se_section,"PRINT%EWALD_INFO",extension=".seLog",error=error) 00705 IF(calculate_forces) THEN 00706 CALL get_qs_env(qs_env=qs_env, force=force, error=error) 00707 00708 ! Allocate atom index for kind 00709 ALLOCATE (atom_of_kind(natoms),STAT=stat) 00710 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00711 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 00712 00713 ! Allocate and zeroing arrays 00714 ALLOCATE ( forces_g(3, forces_g_size), STAT=stat ) 00715 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00716 ALLOCATE ( forces_r(3, natoms), STAT=stat ) 00717 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00718 forces_g = 0.0_dp 00719 forces_r = 0.0_dp 00720 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, se_nonbond_env, cell,& 00721 particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task,& 00722 do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=use_virial, do_efield=.TRUE., & 00723 charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole,& 00724 forces_local=forces_g, forces_glob=forces_r, pv_glob=pv_glob, pv_local=pv_local,& 00725 efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, iw=iw,& 00726 do_debug=.TRUE.,error=error) 00727 ! Only SR force have to be summed up.. the one in g-space are already fully local.. 00728 CALL mp_sum(forces_r, para_env%group) 00729 ELSE 00730 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, se_nonbond_env, cell,& 00731 particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task,& 00732 do_correction_bonded=.FALSE., do_forces=.FALSE., do_stress=.FALSE., do_efield=.TRUE.,& 00733 charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole,& 00734 efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2,& 00735 iw=iw, do_debug=.TRUE.,error=error) 00736 END IF 00737 CALL cp_print_key_finished_output(iw,logger,se_section,"PRINT%EWALD_INFO",error=error) 00738 00739 ! Apply correction only when the Integral Scheme is different from Slater 00740 IF ((se_control%integral_screening/=do_se_IS_slater).AND.(.NOT.debug_this_module)) THEN 00741 CALL build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, calculate_forces,& 00742 store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, virial,& 00743 pv_glob, error) 00744 END IF 00745 00746 ! Virial for the long-range part and correction 00747 IF (use_virial) THEN 00748 ! Sum up contribution of pv_glob on each thread and keep only one copy of pv_local 00749 virial%pv_virial = virial%pv_virial + pv_glob 00750 IF (para_env%mepos==para_env%source) THEN 00751 virial%pv_virial = virial%pv_virial + pv_local 00752 END IF 00753 END IF 00754 00755 ! Debug Statements 00756 IF (debug_this_module) THEN 00757 CALL mp_sum(energy_glob,para_env%group) 00758 WRITE(*,*)"TOTAL ENERGY AFTER EWALD:",energy_local+ energy_glob+ e_neut+ e_self,& 00759 energy_local, energy_glob, e_neut, e_self 00760 END IF 00761 00762 ! Modify the KS matrix and possibly compute derivatives 00763 node = 0 00764 DO ikind=1,nkind 00765 atomic_kind => atomic_kind_set(ikind) 00766 CALL get_atomic_kind(atomic_kind=atomic_kind,se_parameter=se_kind_a) 00767 CALL get_se_param(se_kind_a,defined=defined,natorb=natorb_a) 00768 00769 IF (.NOT.defined .OR. natorb_a < 1) CYCLE 00770 00771 nparticle_local = local_particles%n_el(ikind) 00772 DO ilist=1, nparticle_local 00773 node = node + 1 00774 iatom = local_particles%list(ikind)%array(ilist) 00775 DO ispin = 1, nspins 00776 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(ispin)%matrix,& 00777 row=iatom,col=iatom,BLOCK=ksa_block_a,found=found) 00778 CPPostcondition(ASSOCIATED(ksa_block_a),cp_failure_level,routineP,error,failure) 00779 00780 ! Modify Hamiltonian Matrix accordingly potential, field and electric field gradient 00781 size_1c_int = SIZE(se_kind_a%w_mpole) 00782 DO jint = 1, size_1c_int 00783 tmp = 0.0_dp 00784 mpole => se_kind_a%w_mpole(jint)%mpole 00785 indi = se_orbital_pointer(mpole%indi) 00786 indj = se_orbital_pointer(mpole%indj) 00787 00788 ! Charge 00789 IF (mpole%task(1).AND.task(1)) THEN 00790 tmp = tmp + mpole%c*se_nddo_mpole%efield0(iatom) 00791 END IF 00792 00793 ! Dipole 00794 IF (mpole%task(2).AND.task(2)) THEN 00795 tmp = tmp - DOT_PRODUCT(mpole%d,se_nddo_mpole%efield1(:,iatom)) 00796 END IF 00797 00798 ! Quadrupole 00799 IF (mpole%task(3).AND.task(3)) THEN 00800 tmp = tmp - (1.0_dp/3.0_dp) * SUM(mpole%qc * RESHAPE(se_nddo_mpole%efield2(:,iatom),(/3,3/))) 00801 END IF 00802 00803 ksa_block_a(indi, indj) = ksa_block_a(indi, indj) + tmp 00804 ksa_block_a(indj, indi) = ksa_block_a(indi, indj) 00805 END DO 00806 END DO 00807 00808 ! Nuclear term and forces 00809 IF (task(1)) enuc = enuc + se_kind_a%zeff * se_nddo_mpole%efield0(iatom) 00810 IF (atener) THEN 00811 qs_env%atprop%atecoul(iatom) = qs_env%atprop%atecoul(iatom) +& 00812 0.5_dp*se_kind_a%zeff * se_nddo_mpole%efield0(iatom) 00813 END IF 00814 IF(calculate_forces) THEN 00815 atom_a = atom_of_kind(iatom) 00816 force_a = forces_r(1:3,iatom) + forces_g(1:3,node) 00817 ! Derivatives of the periodic Coulomb Terms 00818 force(ikind)%all_potential(:,atom_a) = force(ikind)%all_potential(:,atom_a) - force_a(:) 00819 END IF 00820 END DO 00821 END DO 00822 ! Sum nuclear energy contribution 00823 CALL mp_sum(enuc, para_env%group) 00824 energy%core_overlap = energy%core_overlap + energy%core_overlap0 + 0.5_dp*enuc 00825 00826 ! Debug Statements 00827 IF (debug_this_module) THEN 00828 WRITE(*,*)"ENUC: ",enuc*0.5_dp 00829 END IF 00830 00831 DO ispin = 1, nspins 00832 CALL cp_dbcsr_sum_replicated( diagmat_ks(ispin)%matrix, error=error ) 00833 CALL cp_dbcsr_distribute(diagmat_ks(ispin)%matrix, error=error) 00834 CALL cp_dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix,& 00835 1.0_dp,1.0_dp,error=error) 00836 END DO 00837 CALL cp_dbcsr_deallocate_matrix_set ( diagmat_p, error=error ) 00838 CALL cp_dbcsr_deallocate_matrix_set ( diagmat_ks, error=error ) 00839 00840 ! Set the Fock matrix contribution to SCP 00841 IF (calculate_forces) THEN 00842 DEALLOCATE(atom_of_kind,stat=stat) 00843 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00844 DEALLOCATE(forces_g,stat=stat) 00845 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00846 DEALLOCATE(forces_r,stat=stat) 00847 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00848 END IF 00849 00850 END IF 00851 00852 CALL timestop(handle) 00853 END SUBROUTINE build_fock_matrix_coulomb_lr 00854 00855 ! ***************************************************************************** 00861 SUBROUTINE build_fock_matrix_coul_lrc (qs_env,ks_matrix,matrix_p,energy,& 00862 calculate_forces, store_int_env,se_nddo_mpole,task,diagmat_p,diagmat_ks,& 00863 virial, pv_glob, error) 00864 00865 TYPE(qs_environment_type), POINTER :: qs_env 00866 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00867 POINTER :: ks_matrix, matrix_p 00868 TYPE(qs_energy_type), POINTER :: energy 00869 LOGICAL, INTENT(in) :: calculate_forces 00870 TYPE(semi_empirical_si_type), POINTER :: store_int_env 00871 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 00872 LOGICAL, DIMENSION(3), INTENT(IN) :: task 00873 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00874 POINTER :: diagmat_p, diagmat_ks 00875 TYPE(virial_type), POINTER :: virial 00876 REAL(KIND=dp), DIMENSION(3, 3), 00877 INTENT(INOUT) :: pv_glob 00878 TYPE(cp_error_type), INTENT(inout) :: error 00879 00880 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coul_lrc', 00881 routineP = moduleN//':'//routineN 00882 00883 INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, itype, jatom, 00884 jkind, natom, natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nspins, 00885 size1, size2, stat 00886 INTEGER, ALLOCATABLE, DIMENSION(:) :: se_natorb 00887 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 00888 LOGICAL :: anag, atener, check, defined, 00889 failure, found, switch, 00890 use_virial 00891 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 00892 REAL(KIND=dp) :: delta, dr1, ecore2, enuc, enuclear, ptens11, ptens12, 00893 ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33 00894 REAL(KIND=dp), DIMENSION(2) :: ecab 00895 REAL(KIND=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b 00896 REAL(KIND=dp), DIMENSION(3) :: force_ab, force_ab0, rij 00897 REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b 00898 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 00899 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2, ksa_block_a, 00900 ksa_block_b, ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, 00901 pb_block_a, pb_block_b 00902 TYPE(atomic_kind_type), DIMENSION(:), 00903 POINTER :: atomic_kind_set 00904 TYPE(atomic_kind_type), POINTER :: atomic_kind 00905 TYPE(cell_type), POINTER :: cell 00906 TYPE(cp_para_env_type), POINTER :: para_env 00907 TYPE(dft_control_type), POINTER :: dft_control 00908 TYPE(neighbor_list_iterator_p_type), 00909 DIMENSION(:), POINTER :: nl_iterator 00910 TYPE(neighbor_list_set_p_type), 00911 DIMENSION(:), POINTER :: sab_lrc 00912 TYPE(particle_type), DIMENSION(:), 00913 POINTER :: particle_set 00914 TYPE(qs_force_type), DIMENSION(:), 00915 POINTER :: force 00916 TYPE(se_int_control_type) :: se_int_control 00917 TYPE(se_taper_type), POINTER :: se_taper 00918 TYPE(semi_empirical_control_type), 00919 POINTER :: se_control 00920 TYPE(semi_empirical_p_type), 00921 DIMENSION(:), POINTER :: se_kind_list 00922 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 00923 00924 failure=.FALSE. 00925 CALL timeset(routineN,handle) 00926 IF ( .NOT. failure ) THEN 00927 NULLIFY(dft_control, cell, force, particle_set, se_control, se_taper,& 00928 efield0, efield1, efield2) 00929 00930 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper,& 00931 para_env=para_env, sab_lrc=sab_lrc, atomic_kind_set=atomic_kind_set, & 00932 particle_set=particle_set, error=error) 00933 00934 ! Parameters 00935 CALL initialize_se_taper(se_taper,lr_corr=.TRUE.,error=error) 00936 se_control => dft_control%qs_control%se_control 00937 anag = se_control%analytical_gradients 00938 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer).AND.calculate_forces 00939 atener = qs_env%atprop%energy 00940 00941 CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3,& 00942 do_ewald_gks=.FALSE., integral_screening=se_control%integral_screening,& 00943 shortrange=.FALSE., max_multipole=se_control%max_multipole,& 00944 pc_coulomb_int=.TRUE.) 00945 00946 nspins=dft_control%nspins 00947 CPPrecondition(ASSOCIATED(matrix_p),cp_failure_level,routineP,error,failure) 00948 CPPrecondition(SIZE(ks_matrix)>0,cp_failure_level,routineP,error,failure) 00949 CPPrecondition(ASSOCIATED(diagmat_p),cp_failure_level,routineP,error,failure) 00950 CPPrecondition(ASSOCIATED(diagmat_ks),cp_failure_level,routineP,error,failure) 00951 00952 nkind = SIZE(atomic_kind_set) 00953 IF(calculate_forces) THEN 00954 CALL get_qs_env(qs_env=qs_env, force=force, error=error) 00955 natom = SIZE (particle_set) 00956 delta = se_control%delta 00957 ! Allocate atom index for kind 00958 ALLOCATE (atom_of_kind(natom),STAT=stat) 00959 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00960 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,atom_of_kind=atom_of_kind) 00961 END IF 00962 00963 ! Allocate arrays for storing partial information on potential, field, field gradient 00964 size1 = SIZE(se_nddo_mpole%efield0) 00965 ALLOCATE (efield0(size1), stat=stat) 00966 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00967 efield0 = 0.0_dp 00968 size1 = SIZE(se_nddo_mpole%efield1,1) 00969 size2 = SIZE(se_nddo_mpole%efield1,2) 00970 ALLOCATE (efield1(size1,size2), stat=stat) 00971 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00972 efield1 = 0.0_dp 00973 size1 = SIZE(se_nddo_mpole%efield2,1) 00974 size2 = SIZE(se_nddo_mpole%efield2,2) 00975 ALLOCATE (efield2(size1,size2), stat=stat) 00976 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00977 efield2 = 0.0_dp 00978 00979 ! Initialize if virial is requested 00980 IF (use_virial) THEN 00981 ptens11 = 0.0_dp ; ptens12 = 0.0_dp ; ptens13 = 0.0_dp 00982 ptens21 = 0.0_dp ; ptens22 = 0.0_dp ; ptens23 = 0.0_dp 00983 ptens31 = 0.0_dp ; ptens32 = 0.0_dp ; ptens33 = 0.0_dp 00984 END IF 00985 00986 ! Start of the loop for the correction of the pair interactions 00987 ecore2 = 0.0_dp 00988 enuclear = 0.0_dp 00989 itype = get_se_type(dft_control%qs_control%method_id) 00990 00991 ALLOCATE (se_defined(nkind),se_kind_list(nkind),se_natorb(nkind),STAT=stat) 00992 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00993 DO ikind=1,nkind 00994 atomic_kind => atomic_kind_set(ikind) 00995 CALL get_atomic_kind(atomic_kind=atomic_kind,se_parameter=se_kind_a) 00996 se_kind_list(ikind)%se_param => se_kind_a 00997 CALL get_se_param(se_kind_a,defined=defined,natorb=natorb_a) 00998 se_defined(ikind) = (defined .AND. natorb_a >= 1) 00999 se_natorb(ikind) = natorb_a 01000 END DO 01001 01002 CALL neighbor_list_iterator_create(nl_iterator,sab_lrc) 01003 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 01004 CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,iatom=iatom,jatom=jatom,inode=inode,r=rij) 01005 IF (.NOT.se_defined(ikind)) CYCLE 01006 IF (.NOT.se_defined(jkind)) CYCLE 01007 se_kind_a => se_kind_list(ikind)%se_param 01008 se_kind_b => se_kind_list(jkind)%se_param 01009 natorb_a = se_natorb(ikind) 01010 natorb_b = se_natorb(jkind) 01011 natorb_a2 = natorb_a**2 01012 natorb_b2 = natorb_b**2 01013 01014 IF (inode==1) THEN 01015 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(1)%matrix,& 01016 row=iatom,col=iatom,BLOCK=pa_block_a,found=found) 01017 CPPostcondition(ASSOCIATED(pa_block_a),cp_failure_level,routineP,error,failure) 01018 check = (SIZE(pa_block_a,1)==natorb_a).AND.(SIZE(pa_block_a,2)==natorb_a) 01019 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01020 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix,& 01021 row=iatom,col=iatom,BLOCK=ksa_block_a,found=found) 01022 CPPostcondition(ASSOCIATED(ksa_block_a),cp_failure_level,routineP,error,failure) 01023 p_block_tot_a(1:natorb_a,1:natorb_a) = 2.0_dp * pa_block_a 01024 pa_a(1:natorb_a2) = RESHAPE(pa_block_a,(/natorb_a2/)) 01025 IF ( nspins == 2 ) THEN 01026 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(2)%matrix,& 01027 row=iatom,col=iatom,BLOCK=pa_block_b,found=found) 01028 CPPostcondition(ASSOCIATED(pa_block_b),cp_failure_level,routineP,error,failure) 01029 check = (SIZE(pa_block_b,1)==natorb_a).AND.(SIZE(pa_block_b,2)==natorb_a) 01030 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01031 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix,& 01032 row=iatom,col=iatom,BLOCK=ksa_block_b,found=found) 01033 CPPostcondition(ASSOCIATED(ksa_block_b),cp_failure_level,routineP,error,failure) 01034 p_block_tot_a(1:natorb_a,1:natorb_a) = pa_block_a + pa_block_b 01035 pa_b(1:natorb_a2) = RESHAPE(pa_block_b,(/natorb_a2/)) 01036 END IF 01037 01038 ENDIF 01039 01040 dr1 = DOT_PRODUCT(rij,rij) 01041 IF ( dr1 > rij_threshold ) THEN 01042 ! Determine the order of the atoms, and in case switch them.. 01043 IF (iatom <= jatom) THEN 01044 switch = .FALSE. 01045 ELSE 01046 switch = .TRUE. 01047 END IF 01048 01049 ! Point-like interaction corrections 01050 CALL se_coulomb_ij_interaction (iatom, jatom, task, do_forces=calculate_forces,& 01051 do_efield=.TRUE., do_stress=use_virial, charges=se_nddo_mpole%charge, & 01052 dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, & 01053 force_ab=force_ab0, efield0=efield0, efield1=efield1, efield2=efield2,& 01054 rab2=dr1, rab=rij, ptens11=ptens11, ptens12=ptens12, ptens13=ptens13,& 01055 ptens21=ptens21, ptens22=ptens22, ptens23=ptens23, ptens31=ptens31,& 01056 ptens32=ptens32, ptens33=ptens33, error=error) 01057 01058 ! Retrieve blocks for KS and P 01059 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(1)%matrix,& 01060 row=jatom,col=jatom,BLOCK=pb_block_a,found=found) 01061 CPPostcondition(ASSOCIATED(pb_block_a),cp_failure_level,routineP,error,failure) 01062 check = (SIZE(pb_block_a,1)==natorb_b).AND.(SIZE(pb_block_a,2)==natorb_b) 01063 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01064 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix,& 01065 row=jatom,col=jatom,BLOCK=ksb_block_a,found=found) 01066 CPPostcondition(ASSOCIATED(ksb_block_a),cp_failure_level,routineP,error,failure) 01067 p_block_tot_b(1:natorb_b,1:natorb_b) = 2.0_dp * pb_block_a 01068 pb_a(1:natorb_b2) = RESHAPE(pb_block_a,(/natorb_b2/)) 01069 ! Handle more than one configuration 01070 IF ( nspins == 2 ) THEN 01071 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(2)%matrix,& 01072 row=jatom,col=jatom,BLOCK=pb_block_b,found=found) 01073 CPPostcondition(ASSOCIATED(pb_block_b),cp_failure_level,routineP,error,failure) 01074 check = (SIZE(pb_block_b,1)==natorb_b).AND.(SIZE(pb_block_b,2)==natorb_b) 01075 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01076 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix,& 01077 row=jatom,col=jatom,BLOCK=ksb_block_b,found=found) 01078 CPPostcondition(ASSOCIATED(ksb_block_b),cp_failure_level,routineP,error,failure) 01079 check = (SIZE(pb_block_a,1)==SIZE(pb_block_b,1)).AND.(SIZE(pb_block_a,2)==SIZE(pb_block_b,2)) 01080 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01081 p_block_tot_b(1:natorb_b,1:natorb_b) = pb_block_a + pb_block_b 01082 pb_b(1:natorb_b2) = RESHAPE(pb_block_b,(/natorb_b2/)) 01083 END IF 01084 01085 SELECT CASE (dft_control%qs_control%method_id) 01086 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pdg,& 01087 do_method_rm1, do_method_mndod) 01088 ! Evaluate nuclear contribution.. 01089 CALL corecore_el (se_kind_a,se_kind_b,rij,enuc=enuc,itype=itype,anag=anag,& 01090 se_int_control=se_int_control, se_taper=se_taper, error=error) 01091 enuclear = enuclear + enuc 01092 01093 ! Two-centers One-electron terms 01094 IF ( nspins == 1 ) THEN 01095 ecab = 0._dp 01096 CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_a, ksb_block_a,& 01097 pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control,& 01098 se_taper=se_taper, store_int_env=store_int_env, error=error) 01099 ecore2 = ecore2 + ecab(1) + ecab(2) 01100 ELSE IF ( nspins == 2 ) THEN 01101 ecab = 0._dp 01102 CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_a, ksb_block_a,& 01103 pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag,& 01104 se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env,& 01105 error=error) 01106 CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_b, ksb_block_b,& 01107 pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control,& 01108 se_taper=se_taper, store_int_env=store_int_env, error=error) 01109 ecore2 = ecore2 + ecab(1) + ecab(2) 01110 END IF 01111 IF (atener) THEN 01112 qs_env%atprop%atecoul(iatom) = qs_env%atprop%atecoul(iatom) + ecab(1) + 0.5_dp*enuc 01113 qs_env%atprop%atecoul(jatom) = qs_env%atprop%atecoul(jatom) + ecab(2) + 0.5_dp*enuc 01114 END IF 01115 ! Coulomb Terms 01116 IF ( nspins == 1 ) THEN 01117 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,& 01118 ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 01119 store_int_env=store_int_env, error=error) 01120 ELSE IF ( nspins == 2 ) THEN 01121 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,& 01122 ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 01123 store_int_env=store_int_env, error=error) 01124 01125 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b,& 01126 ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,& 01127 store_int_env=store_int_env, error=error) 01128 END IF 01129 01130 IF(calculate_forces) THEN 01131 atom_a = atom_of_kind(iatom) 01132 atom_b = atom_of_kind(jatom) 01133 01134 ! Evaluate nuclear contribution.. 01135 CALL dcorecore_el (se_kind_a,se_kind_b,rij,denuc=force_ab,itype=itype,delta=delta,& 01136 anag=anag,se_int_control=se_int_control,se_taper=se_taper,error=error) 01137 01138 ! Derivatives of the Two-centre One-electron terms 01139 IF ( nspins == 1 ) THEN 01140 CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_a, pb_a, itype=itype, anag=anag,& 01141 se_int_control=se_int_control, se_taper=se_taper, force=force_ab,& 01142 delta=delta, error=error) 01143 ELSE IF ( nspins == 2 ) THEN 01144 CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_block_a, pb_block_a, itype=itype, anag=anag,& 01145 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 01146 error=error) 01147 CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_b, pb_b, itype=itype, anag=anag,& 01148 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 01149 error=error) 01150 END IF 01151 IF (use_virial) THEN 01152 CALL virial_pair_force ( virial%pv_virial, -1.0_dp, force_ab, rij, error) 01153 END IF 01154 force_ab = force_ab + force_ab0 01155 01156 ! Sum up force components 01157 force(ikind)%all_potential(1,atom_a) = force(ikind)%all_potential(1,atom_a) - force_ab(1) 01158 force(jkind)%all_potential(1,atom_b) = force(jkind)%all_potential(1,atom_b) + force_ab(1) 01159 01160 force(ikind)%all_potential(2,atom_a) = force(ikind)%all_potential(2,atom_a) - force_ab(2) 01161 force(jkind)%all_potential(2,atom_b) = force(jkind)%all_potential(2,atom_b) + force_ab(2) 01162 01163 force(ikind)%all_potential(3,atom_a) = force(ikind)%all_potential(3,atom_a) - force_ab(3) 01164 force(jkind)%all_potential(3,atom_b) = force(jkind)%all_potential(3,atom_b) + force_ab(3) 01165 01166 ! Derivatives of the Coulomb Terms 01167 force_ab = 0._dp 01168 IF ( nspins == 1 ) THEN 01169 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp,& 01170 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 01171 error=error) 01172 ELSE IF ( nspins == 2 ) THEN 01173 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,& 01174 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 01175 error=error) 01176 01177 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,& 01178 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,& 01179 error=error) 01180 END IF 01181 IF ( switch ) THEN 01182 force_ab(1) = -force_ab(1) 01183 force_ab(2) = -force_ab(2) 01184 force_ab(3) = -force_ab(3) 01185 END IF 01186 IF (use_virial) THEN 01187 CALL virial_pair_force ( virial%pv_virial, -1.0_dp, force_ab, rij, error) 01188 END IF 01189 01190 ! Sum up force components 01191 force(ikind)%rho_elec(1,atom_a) = force(ikind)%rho_elec(1,atom_a) - force_ab(1) 01192 force(jkind)%rho_elec(1,atom_b) = force(jkind)%rho_elec(1,atom_b) + force_ab(1) 01193 01194 force(ikind)%rho_elec(2,atom_a) = force(ikind)%rho_elec(2,atom_a) - force_ab(2) 01195 force(jkind)%rho_elec(2,atom_b) = force(jkind)%rho_elec(2,atom_b) + force_ab(2) 01196 01197 force(ikind)%rho_elec(3,atom_a) = force(ikind)%rho_elec(3,atom_a) - force_ab(3) 01198 force(jkind)%rho_elec(3,atom_b) = force(jkind)%rho_elec(3,atom_b) + force_ab(3) 01199 END IF 01200 CASE DEFAULT 01201 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01202 END SELECT 01203 END IF 01204 END DO 01205 CALL neighbor_list_iterator_release(nl_iterator) 01206 01207 DEALLOCATE(se_kind_list,se_defined,se_natorb,stat=stat) 01208 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01209 01210 ! Sum-up Virial constribution (long-range correction) 01211 IF (use_virial) THEN 01212 pv_glob(1,1) = pv_glob(1,1) + ptens11 01213 pv_glob(1,2) = pv_glob(1,2) + (ptens12+ptens21)*0.5_dp 01214 pv_glob(1,3) = pv_glob(1,3) + (ptens13+ptens31)*0.5_dp 01215 pv_glob(2,1) = pv_glob(1,2) 01216 pv_glob(2,2) = pv_glob(2,2) + ptens22 01217 pv_glob(2,3) = pv_glob(2,3) + (ptens23+ptens32)*0.5_dp 01218 pv_glob(3,1) = pv_glob(1,3) 01219 pv_glob(3,2) = pv_glob(2,3) 01220 pv_glob(3,3) = pv_glob(3,3) + ptens33 01221 END IF 01222 01223 IF (calculate_forces) THEN 01224 DEALLOCATE(atom_of_kind,stat=stat) 01225 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01226 END IF 01227 01228 ! collect information on potential, field, field gradient 01229 CALL mp_sum(efield0, para_env%group) 01230 CALL mp_sum(efield1, para_env%group) 01231 CALL mp_sum(efield2, para_env%group) 01232 se_nddo_mpole%efield0 = se_nddo_mpole%efield0 - efield0 01233 se_nddo_mpole%efield1 = se_nddo_mpole%efield1 - efield1 01234 se_nddo_mpole%efield2 = se_nddo_mpole%efield2 - efield2 01235 ! deallocate working arrays 01236 DEALLOCATE (efield0, stat=stat) 01237 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01238 DEALLOCATE (efield1, stat=stat) 01239 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01240 DEALLOCATE (efield2, stat=stat) 01241 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01242 01243 ! Corrections for two-centers one-electron terms + nuclear terms 01244 CALL mp_sum(enuclear,para_env%group) 01245 CALL mp_sum(ecore2,para_env%group) 01246 energy%hartree = energy%hartree + ecore2 01247 energy%core_overlap = enuclear 01248 END IF 01249 CALL finalize_se_taper(se_taper,error=error) 01250 CALL timestop(handle) 01251 END SUBROUTINE build_fock_matrix_coul_lrc 01252 01253 ! ***************************************************************************** 01260 SUBROUTINE build_fock_matrix_coul_lr_r3(qs_env,ks_matrix,matrix_p,energy,calculate_forces,error) 01261 TYPE(qs_environment_type), POINTER :: qs_env 01262 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01263 POINTER :: ks_matrix, matrix_p 01264 TYPE(qs_energy_type), POINTER :: energy 01265 LOGICAL, INTENT(in) :: calculate_forces 01266 TYPE(cp_error_type), INTENT(inout) :: error 01267 01268 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coul_lr_r3', 01269 routineP = moduleN//':'//routineN 01270 01271 INTEGER :: atom_a, atom_b, ewald_type, handle, iatom, ikind, inode, 01272 ispin, itype, jatom, jkind, natoms, natorb_a, natorb_a2, natorb_b, 01273 natorb_b2, nkind, nlocal_particles, nspins, stat 01274 INTEGER, ALLOCATABLE, DIMENSION(:) :: se_natorb 01275 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 01276 LOGICAL :: anag, atener, check, defined, 01277 failure, found, switch 01278 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 01279 REAL(KIND=dp) :: dr1, ecore2, r2inv, r3inv, 01280 rinv 01281 REAL(KIND=dp), DIMENSION(2) :: ecab 01282 REAL(KIND=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b 01283 REAL(KIND=dp), DIMENSION(3) :: dr3inv, force_ab, rij 01284 REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b 01285 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, 01286 ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, pb_block_a, pb_block_b 01287 TYPE(atomic_kind_type), DIMENSION(:), 01288 POINTER :: atomic_kind_set 01289 TYPE(atomic_kind_type), POINTER :: atomic_kind 01290 TYPE(cell_type), POINTER :: cell 01291 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01292 POINTER :: diagmat_ks, diagmat_p 01293 TYPE(cp_logger_type), POINTER :: logger 01294 TYPE(cp_para_env_type), POINTER :: para_env 01295 TYPE(dft_control_type), POINTER :: dft_control 01296 TYPE(distribution_1d_type), POINTER :: local_particles 01297 TYPE(ewald_environment_type), POINTER :: ewald_env 01298 TYPE(ewald_pw_type), POINTER :: ewald_pw 01299 TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env 01300 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 01301 TYPE(neighbor_list_iterator_p_type), 01302 DIMENSION(:), POINTER :: nl_iterator 01303 TYPE(neighbor_list_set_p_type), 01304 DIMENSION(:), POINTER :: sab_orb 01305 TYPE(particle_type), DIMENSION(:), 01306 POINTER :: particle_set 01307 TYPE(qs_force_type), DIMENSION(:), 01308 POINTER :: force 01309 TYPE(section_vals_type), POINTER :: se_section 01310 TYPE(semi_empirical_control_type), 01311 POINTER :: se_control 01312 TYPE(semi_empirical_p_type), 01313 DIMENSION(:), POINTER :: se_kind_list 01314 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 01315 01316 failure=.FALSE. 01317 CALL timeset(routineN,handle) 01318 01319 CALL get_qs_env(qs_env=qs_env, error=error)!sm->dbcsr 01320 01321 IF ( .NOT. failure ) THEN 01322 NULLIFY(dft_control, cell, force, particle_set, diagmat_ks, & 01323 diagmat_p, local_particles, se_control, ewald_env, ewald_pw, & 01324 se_nddo_mpole, se_nonbond_env, se_section, sab_orb, logger) 01325 01326 logger => cp_error_get_logger(error) 01327 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env,& 01328 atomic_kind_set=atomic_kind_set, particle_set=particle_set, ewald_env=ewald_env,& 01329 local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole,& 01330 se_nonbond_env=se_nonbond_env, sab_orb=sab_orb, error=error) 01331 01332 nlocal_particles = SUM(local_particles%n_el(:)) 01333 natoms = SIZE(particle_set) 01334 CALL ewald_env_get (ewald_env, ewald_type=ewald_type, error=error) 01335 SELECT CASE(ewald_type) 01336 CASE (do_ewald_ewald) 01337 ! Do Nothing 01338 CASE DEFAULT 01339 CALL cp_unimplemented_error(fromWhere=routineP, & 01340 message="Periodic SE implemented only for standard EWALD sums.", & 01341 error=error, error_level=cp_failure_level) 01342 END SELECT 01343 01344 ! Parameters 01345 se_section => section_vals_get_subs_vals(qs_env%input,"DFT%QS%SE",error=error) 01346 se_control => dft_control%qs_control%se_control 01347 anag = se_control%analytical_gradients 01348 atener = qs_env%atprop%energy 01349 01350 nspins=dft_control%nspins 01351 CPPrecondition(ASSOCIATED(matrix_p),cp_failure_level,routineP,error,failure) 01352 CPPrecondition(SIZE(ks_matrix)>0,cp_failure_level,routineP,error,failure) 01353 01354 CALL cp_dbcsr_allocate_matrix_set(diagmat_ks,nspins,error=error) 01355 CALL cp_dbcsr_allocate_matrix_set(diagmat_p,nspins,error=error) 01356 01357 nkind = SIZE(atomic_kind_set) 01358 DO ispin=1,nspins 01359 ! Allocate diagonal block matrices 01360 ALLOCATE(diagmat_p(ispin)%matrix,diagmat_ks(ispin)%matrix)!sm->dbcsr 01361 CALL cp_dbcsr_init(diagmat_p(ispin)%matrix, error=error) 01362 CALL cp_dbcsr_init(diagmat_ks(ispin)%matrix, error=error) 01363 CALL cp_dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix, error=error) 01364 CALL cp_dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, error=error) 01365 CALL cp_dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp, error=error) 01366 CALL cp_dbcsr_replicate_all(diagmat_p(ispin)%matrix, error=error) 01367 CALL cp_dbcsr_replicate_all(diagmat_ks(ispin)%matrix, error=error) 01368 END DO 01369 01370 ! Possibly compute forces 01371 IF(calculate_forces) THEN 01372 CALL get_qs_env(qs_env=qs_env,force=force,error=error) 01373 ALLOCATE (atom_of_kind(natoms),STAT=stat) 01374 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01375 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,atom_of_kind=atom_of_kind) 01376 END IF 01377 itype = get_se_type(dft_control%qs_control%method_id) 01378 01379 ecore2 = 0.0_dp 01380 ! Real space part of the 1/R^3 sum 01381 ALLOCATE (se_defined(nkind),se_kind_list(nkind),se_natorb(nkind),STAT=stat) 01382 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01383 DO ikind=1,nkind 01384 atomic_kind => atomic_kind_set(ikind) 01385 CALL get_atomic_kind(atomic_kind=atomic_kind,se_parameter=se_kind_a) 01386 se_kind_list(ikind)%se_param => se_kind_a 01387 CALL get_se_param(se_kind_a,defined=defined,natorb=natorb_a) 01388 se_defined(ikind) = (defined .AND. natorb_a >= 1) 01389 se_natorb(ikind) = natorb_a 01390 END DO 01391 01392 CALL neighbor_list_iterator_create(nl_iterator,sab_orb) 01393 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 01394 CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,iatom=iatom,jatom=jatom,inode=inode,r=rij) 01395 IF (.NOT.se_defined(ikind)) CYCLE 01396 IF (.NOT.se_defined(jkind)) CYCLE 01397 se_kind_a => se_kind_list(ikind)%se_param 01398 se_kind_b => se_kind_list(jkind)%se_param 01399 natorb_a = se_natorb(ikind) 01400 natorb_b = se_natorb(jkind) 01401 natorb_a2 = natorb_a**2 01402 natorb_b2 = natorb_b**2 01403 01404 IF (inode==1) THEN 01405 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(1)%matrix,& 01406 row=iatom,col=iatom,BLOCK=pa_block_a,found=found) 01407 CPPostcondition(ASSOCIATED(pa_block_a),cp_failure_level,routineP,error,failure) 01408 check = (SIZE(pa_block_a,1)==natorb_a).AND.(SIZE(pa_block_a,2)==natorb_a) 01409 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01410 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix,& 01411 row=iatom,col=iatom,BLOCK=ksa_block_a,found=found) 01412 CPPostcondition(ASSOCIATED(ksa_block_a),cp_failure_level,routineP,error,failure) 01413 p_block_tot_a(1:natorb_a,1:natorb_a) = 2.0_dp * pa_block_a 01414 pa_a(1:natorb_a2) = RESHAPE(pa_block_a,(/natorb_a2/)) 01415 IF ( nspins == 2 ) THEN 01416 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(2)%matrix,& 01417 row=iatom,col=iatom,BLOCK=pa_block_b,found=found) 01418 CPPostcondition(ASSOCIATED(pa_block_b),cp_failure_level,routineP,error,failure) 01419 check = (SIZE(pa_block_b,1)==natorb_a).AND.(SIZE(pa_block_b,2)==natorb_a) 01420 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01421 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix,& 01422 row=iatom,col=iatom,BLOCK=ksa_block_b,found=found) 01423 CPPostcondition(ASSOCIATED(ksa_block_b),cp_failure_level,routineP,error,failure) 01424 p_block_tot_a(1:natorb_a,1:natorb_a) = pa_block_a + pa_block_b 01425 pa_b(1:natorb_a2) = RESHAPE(pa_block_b,(/natorb_a2/)) 01426 END IF 01427 END IF 01428 01429 dr1 = DOT_PRODUCT(rij,rij) 01430 IF ( dr1 > rij_threshold ) THEN 01431 ! Determine the order of the atoms, and in case switch them.. 01432 IF (iatom <= jatom) THEN 01433 switch = .FALSE. 01434 ELSE 01435 switch = .TRUE. 01436 END IF 01437 ! Retrieve blocks for KS and P 01438 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(1)%matrix,& 01439 row=jatom,col=jatom,BLOCK=pb_block_a,found=found) 01440 CPPostcondition(ASSOCIATED(pb_block_a),cp_failure_level,routineP,error,failure) 01441 check = (SIZE(pb_block_a,1)==natorb_b).AND.(SIZE(pb_block_a,2)==natorb_b) 01442 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01443 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix,& 01444 row=jatom,col=jatom,BLOCK=ksb_block_a,found=found) 01445 CPPostcondition(ASSOCIATED(ksb_block_a),cp_failure_level,routineP,error,failure) 01446 p_block_tot_b(1:natorb_b,1:natorb_b) = 2.0_dp * pb_block_a 01447 pb_a(1:natorb_b2) = RESHAPE(pb_block_a,(/natorb_b2/)) 01448 ! Handle more than one configuration 01449 IF ( nspins == 2 ) THEN 01450 CALL cp_dbcsr_get_block_p(matrix=diagmat_p(2)%matrix,& 01451 row=jatom,col=jatom,BLOCK=pb_block_b,found=found) 01452 CPPostcondition(ASSOCIATED(pb_block_b),cp_failure_level,routineP,error,failure) 01453 check = (SIZE(pb_block_b,1)==natorb_b).AND.(SIZE(pb_block_b,2)==natorb_b) 01454 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01455 CALL cp_dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix,& 01456 row=jatom,col=jatom,BLOCK=ksb_block_b,found=found) 01457 CPPostcondition(ASSOCIATED(ksb_block_b),cp_failure_level,routineP,error,failure) 01458 check = (SIZE(pb_block_a,1)==SIZE(pb_block_b,1)).AND.(SIZE(pb_block_a,2)==SIZE(pb_block_b,2)) 01459 CPPostcondition(check,cp_failure_level,routineP,error,failure) 01460 p_block_tot_b(1:natorb_b,1:natorb_b) = pb_block_a + pb_block_b 01461 pb_b(1:natorb_b2) = RESHAPE(pb_block_b,(/natorb_b2/)) 01462 END IF 01463 01464 SELECT CASE (dft_control%qs_control%method_id) 01465 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pdg,& 01466 do_method_rm1, do_method_mndod, do_method_pnnl) 01467 01468 ! Pre-compute some quantities.. 01469 r2inv = 1.0_dp/dr1 01470 rinv = SQRT(r2inv) 01471 r3inv = rinv**3 01472 01473 ! Two-centers One-electron terms 01474 IF ( nspins == 1 ) THEN 01475 ecab = 0._dp 01476 CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_a, pb_a,& 01477 ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b,& 01478 e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv, error=error) 01479 ecore2 = ecore2 + ecab(1) + ecab(2) 01480 ELSE IF ( nspins == 2 ) THEN 01481 ecab = 0._dp 01482 CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_block_a,& 01483 pb_block_a, ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b,& 01484 e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv, error=error) 01485 01486 CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_b, ksb_block_b, pa_b, pb_b,& 01487 ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b,& 01488 e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv, error=error) 01489 ecore2 = ecore2 + ecab(1) + ecab(2) 01490 END IF 01491 IF (atener) THEN 01492 qs_env%atprop%atecoul(iatom) = qs_env%atprop%atecoul(iatom) + ecab(1) 01493 qs_env%atprop%atecoul(jatom) = qs_env%atprop%atecoul(jatom) + ecab(2) 01494 END IF 01495 ! Coulomb Terms 01496 IF ( nspins == 1 ) THEN 01497 CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,& 01498 ksb_block_a, factor=0.5_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv,& 01499 error=error) 01500 ELSE IF ( nspins == 2 ) THEN 01501 CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,& 01502 ksb_block_a, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv,& 01503 error=error) 01504 01505 CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_b, p_block_tot_b,& 01506 ksb_block_b, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv,& 01507 error=error) 01508 END IF 01509 01510 ! Compute forces if requested 01511 IF(calculate_forces) THEN 01512 dr3inv = -3.0_dp*rij*r3inv*r2inv 01513 atom_a = atom_of_kind(iatom) 01514 atom_b = atom_of_kind(jatom) 01515 01516 force_ab = 0.0_dp 01517 ! Derivatives of the One-centre One-electron terms 01518 IF ( nspins == 1 ) THEN 01519 CALL dfock2_1el_r3(se_kind_a,se_kind_b, dr3inv, pa_a, pb_a, force_ab,& 01520 se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a,& 01521 error=error) 01522 ELSE IF ( nspins == 2 ) THEN 01523 CALL dfock2_1el_r3(se_kind_a,se_kind_b, dr3inv, pa_block_a, pb_block_a, force_ab,& 01524 se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a,& 01525 error=error) 01526 01527 CALL dfock2_1el_r3(se_kind_a,se_kind_b, dr3inv, pa_b, pb_b, force_ab,& 01528 se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a,& 01529 error=error) 01530 END IF 01531 01532 ! Sum up force components 01533 force(ikind)%all_potential(1,atom_a) = force(ikind)%all_potential(1,atom_a) - force_ab(1) 01534 force(jkind)%all_potential(1,atom_b) = force(jkind)%all_potential(1,atom_b) + force_ab(1) 01535 01536 force(ikind)%all_potential(2,atom_a) = force(ikind)%all_potential(2,atom_a) - force_ab(2) 01537 force(jkind)%all_potential(2,atom_b) = force(jkind)%all_potential(2,atom_b) + force_ab(2) 01538 01539 force(ikind)%all_potential(3,atom_a) = force(ikind)%all_potential(3,atom_a) - force_ab(3) 01540 force(jkind)%all_potential(3,atom_b) = force(jkind)%all_potential(3,atom_b) + force_ab(3) 01541 01542 ! Derivatives of the Coulomb Terms 01543 force_ab = 0.0_dp 01544 IF ( nspins == 1 ) THEN 01545 CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp,& 01546 w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab, error=error) 01547 ELSE IF ( nspins == 2 ) THEN 01548 CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,& 01549 w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab, error=error) 01550 01551 CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,& 01552 w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab, error=error) 01553 END IF 01554 01555 ! Sum up force components 01556 force(ikind)%rho_elec(1,atom_a) = force(ikind)%rho_elec(1,atom_a) - force_ab(1) 01557 force(jkind)%rho_elec(1,atom_b) = force(jkind)%rho_elec(1,atom_b) + force_ab(1) 01558 01559 force(ikind)%rho_elec(2,atom_a) = force(ikind)%rho_elec(2,atom_a) - force_ab(2) 01560 force(jkind)%rho_elec(2,atom_b) = force(jkind)%rho_elec(2,atom_b) + force_ab(2) 01561 01562 force(ikind)%rho_elec(3,atom_a) = force(ikind)%rho_elec(3,atom_a) - force_ab(3) 01563 force(jkind)%rho_elec(3,atom_b) = force(jkind)%rho_elec(3,atom_b) + force_ab(3) 01564 END IF 01565 CASE DEFAULT 01566 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01567 END SELECT 01568 END IF 01569 END DO 01570 CALL neighbor_list_iterator_release(nl_iterator) 01571 01572 DEALLOCATE(se_kind_list,se_defined,se_natorb,stat=stat) 01573 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01574 01575 DO ispin = 1, nspins 01576 CALL cp_dbcsr_sum_replicated( diagmat_ks(ispin)%matrix, error=error ) 01577 CALL cp_dbcsr_distribute(diagmat_ks(ispin)%matrix, error=error) 01578 CALL cp_dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix,& 01579 1.0_dp,1.0_dp,error=error) 01580 END DO 01581 CALL cp_dbcsr_deallocate_matrix_set ( diagmat_p, error=error ) 01582 CALL cp_dbcsr_deallocate_matrix_set ( diagmat_ks, error=error ) 01583 01584 IF (calculate_forces) THEN 01585 DEALLOCATE(atom_of_kind,stat=stat) 01586 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01587 END IF 01588 01589 ! Two-centers one-electron terms 01590 CALL mp_sum(ecore2,para_env%group) 01591 energy%hartree = energy%hartree + ecore2 01592 END IF 01593 01594 CALL timestop(handle) 01595 END SUBROUTINE build_fock_matrix_coul_lr_r3 01596 01597 END MODULE se_fock_matrix_coulomb_mpi 01598
1.7.3