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