|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE ewalds_multipole 00013 USE atomic_kind_types, ONLY: atomic_kind_type 00014 USE bibliography, ONLY: Aguado2003,& 00015 Laino2008,& 00016 cite_reference 00017 USE cell_types, ONLY: cell_type,& 00018 pbc 00019 USE damping_dipole_types, ONLY: damping_type,& 00020 no_damping,& 00021 tang_toennies 00022 USE dg_rho0_types, ONLY: dg_rho0_type 00023 USE dg_types, ONLY: dg_get,& 00024 dg_type 00025 USE distribution_1d_types, ONLY: distribution_1d_type 00026 USE erf_fn, ONLY: erf,& 00027 erfc 00028 USE ewald_environment_types, ONLY: ewald_env_get,& 00029 ewald_environment_type 00030 USE ewald_pw_types, ONLY: ewald_pw_get,& 00031 ewald_pw_type 00032 USE f77_blas 00033 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 00034 neighbor_kind_pairs_type 00035 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 00036 fist_nonbond_env_type,& 00037 pos_type 00038 USE input_section_types, ONLY: section_vals_type 00039 USE kinds, ONLY: dp 00040 USE mathconstants, ONLY: fourpi,& 00041 oorootpi,& 00042 pi,& 00043 sqrthalf 00044 USE mathlib, ONLY: matvec_3x3 00045 USE message_passing, ONLY: mp_sum 00046 USE particle_types, ONLY: particle_type 00047 USE pw_grid_types, ONLY: pw_grid_type 00048 USE pw_pool_types, ONLY: pw_pool_type 00049 USE structure_factor_types, ONLY: structure_factor_type 00050 USE structure_factors, ONLY: structure_factor_allocate,& 00051 structure_factor_deallocate,& 00052 structure_factor_evaluate 00053 USE timings, ONLY: timeset,& 00054 timestop 00055 #include "cp_common_uses.h" 00056 00057 IMPLICIT NONE 00058 PRIVATE 00059 #include "ewalds_multipole_debug.h" 00060 00061 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. 00062 LOGICAL, PRIVATE, PARAMETER :: debug_r_space =.FALSE. 00063 LOGICAL, PRIVATE, PARAMETER :: debug_g_space =.FALSE. 00064 LOGICAL, PRIVATE, PARAMETER :: debug_e_field =.FALSE. 00065 LOGICAL, PRIVATE, PARAMETER :: debug_e_field_en =.FALSE. 00066 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewalds_multipole' 00067 00068 PUBLIC :: ewald_multipole_evaluate 00069 00070 CONTAINS 00071 00072 ! ***************************************************************************** 00086 RECURSIVE SUBROUTINE ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env,& 00087 cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self,& 00088 task, do_correction_bonded, do_forces, do_stress, do_efield, radii, charges, dipoles,& 00089 quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1,& 00090 efield2, iw, do_debug, atomic_kind_set, mm_section, error ) 00091 TYPE(ewald_environment_type), POINTER :: ewald_env 00092 TYPE(ewald_pw_type), POINTER :: ewald_pw 00093 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00094 TYPE(cell_type), POINTER :: cell 00095 TYPE(particle_type), POINTER :: particle_set(:) 00096 TYPE(distribution_1d_type), POINTER :: local_particles 00097 REAL(KIND=dp), INTENT(INOUT) :: energy_local, energy_glob 00098 REAL(KIND=dp), INTENT(OUT) :: e_neut, e_self 00099 LOGICAL, DIMENSION(3), INTENT(IN) :: task 00100 LOGICAL, INTENT(IN) :: do_correction_bonded, 00101 do_forces, do_stress, 00102 do_efield 00103 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00104 POINTER :: radii, charges 00105 REAL(KIND=dp), DIMENSION(:, :), 00106 OPTIONAL, POINTER :: dipoles 00107 REAL(KIND=dp), DIMENSION(:, :, :), 00108 OPTIONAL, POINTER :: quadrupoles 00109 REAL(KIND=dp), DIMENSION(:, :), 00110 INTENT(INOUT), OPTIONAL :: forces_local, forces_glob, 00111 pv_local, pv_glob 00112 REAL(KIND=dp), DIMENSION(:), 00113 INTENT(OUT), OPTIONAL :: efield0 00114 REAL(KIND=dp), DIMENSION(:, :), 00115 INTENT(OUT), OPTIONAL :: efield1, efield2 00116 INTEGER, INTENT(IN) :: iw 00117 LOGICAL, INTENT(IN) :: do_debug 00118 TYPE(atomic_kind_type), DIMENSION(:), 00119 OPTIONAL, POINTER :: atomic_kind_set 00120 TYPE(section_vals_type), OPTIONAL, 00121 POINTER :: mm_section 00122 TYPE(cp_error_type), INTENT(inout) :: error 00123 00124 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_evaluate', 00125 routineP = moduleN//':'//routineN 00126 00127 INTEGER :: group, handle, i, j, size1, 00128 size2, stat 00129 LOGICAL :: check_debug, check_efield, 00130 check_forces, do_task(3), 00131 failure 00132 LOGICAL, DIMENSION(3, 3) :: my_task 00133 REAL(KIND=dp) :: e_bonded, e_bonded_t, 00134 e_rspace, e_rspace_t, 00135 energy_glob_t 00136 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0_lr, efield0_sr 00137 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1_lr, efield1_sr, 00138 efield2_lr, efield2_sr 00139 00140 failure = .FALSE. 00141 CALL cite_reference(Aguado2003) 00142 CALL cite_reference(Laino2008) 00143 CALL timeset(routineN,handle) 00144 CPPostcondition(ASSOCIATED(nonbond_env),cp_failure_level,routineP,error,failure) 00145 check_debug = (debug_this_module.OR.debug_r_space.OR.debug_g_space.OR.debug_e_field.OR.debug_e_field_en)& 00146 .EQV.debug_this_module 00147 CPPostcondition(check_debug,cp_failure_level,routineP,error,failure) 00148 check_forces = do_forces.EQV.(PRESENT(forces_local).AND.PRESENT(forces_glob)) 00149 CPPostcondition(check_forces,cp_failure_level,routineP,error,failure) 00150 check_efield = do_efield.EQV.(PRESENT(efield0).OR.PRESENT(efield1).OR.PRESENT(efield2)) 00151 CPPostcondition(check_efield,cp_failure_level,routineP,error,failure) 00152 ! Debugging this module 00153 IF (debug_this_module.AND.do_debug) THEN 00154 ! Debug specifically real space part 00155 IF (debug_r_space) THEN 00156 CALL debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, & 00157 particle_set, local_particles, iw, debug_r_space, error) 00158 STOP "Debug Multipole Requested: Real Part!" 00159 END IF 00160 ! Debug electric fields and gradients as pure derivatives 00161 IF (debug_e_field) THEN 00162 CPPostcondition(PRESENT(atomic_kind_set),cp_failure_level,routineP,error,failure) 00163 CPPostcondition(PRESENT(mm_section),cp_failure_level,routineP,error,failure) 00164 CALL debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env,& 00165 cell, particle_set, local_particles, radii, charges, dipoles,& 00166 quadrupoles, task, iw, atomic_kind_set, mm_section, error) 00167 STOP "Debug Multipole Requested: POT+EFIELDS+GRAD!" 00168 END IF 00169 ! Debug the potential, electric fields and electric fields gradient in oder 00170 ! to retrieve the correct energy 00171 IF (debug_e_field_en) THEN 00172 CALL debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env,& 00173 cell, particle_set, local_particles, radii, charges, dipoles,& 00174 quadrupoles, task, iw, error) 00175 STOP "Debug Multipole Requested: POT+EFIELDS+GRAD to give the correct energy!!" 00176 END IF 00177 END IF 00178 00179 ! Setup the tasks (needed to skip useless parts in the real-space term) 00180 do_task = task 00181 DO i = 1, 3 00182 IF (do_task(i)) THEN 00183 SELECT CASE(i) 00184 CASE(1) 00185 do_task(1) = ANY(charges/=0.0_dp) 00186 CASE(2) 00187 do_task(2) = ANY(dipoles/=0.0_dp) 00188 CASE(3) 00189 do_task(3) = ANY(quadrupoles/=0.0_dp) 00190 END SELECT 00191 END IF 00192 END DO 00193 DO i = 1,3 00194 DO j =i,3 00195 my_task(j,i) = do_task(i).AND.do_task(j) 00196 my_task(i,j) = my_task(j,i) 00197 END DO 00198 END DO 00199 00200 ! Allocate arrays for the evaluation of the potential, fields and electrostatic field gradients 00201 NULLIFY(efield0_sr, efield0_lr, efield1_sr, efield1_lr, efield2_sr, efield2_lr) 00202 IF (do_efield) THEN 00203 IF (PRESENT(efield0)) THEN 00204 size1 = SIZE(efield0) 00205 ALLOCATE (efield0_sr(size1), stat=stat) 00206 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00207 ALLOCATE (efield0_lr(size1), stat=stat) 00208 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00209 efield0_sr = 0.0_dp 00210 efield0_lr = 0.0_dp 00211 END IF 00212 IF (PRESENT(efield1)) THEN 00213 size1 = SIZE(efield1,1) 00214 size2 = SIZE(efield1,2) 00215 ALLOCATE (efield1_sr(size1,size2), stat=stat) 00216 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00217 ALLOCATE (efield1_lr(size1,size2), stat=stat) 00218 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00219 efield1_sr = 0.0_dp 00220 efield1_lr = 0.0_dp 00221 END IF 00222 IF (PRESENT(efield2)) THEN 00223 size1 = SIZE(efield2,1) 00224 size2 = SIZE(efield2,2) 00225 ALLOCATE (efield2_sr(size1,size2), stat=stat) 00226 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00227 ALLOCATE (efield2_lr(size1,size2), stat=stat) 00228 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00229 efield2_sr = 0.0_dp 00230 efield2_lr = 0.0_dp 00231 END IF 00232 END IF 00233 00234 e_rspace = 0.0_dp 00235 e_bonded = 0.0_dp 00236 IF ((.NOT.debug_g_space) .AND. (nonbond_env%do_nonbonded)) THEN 00237 ! Compute the Real Space (Short-Range) part of the Ewald sum. 00238 ! This contribution is only added when the nonbonded flag in the input 00239 ! is set, because these contributions depend. the neighborlists. 00240 CALL ewald_multipole_SR (nonbond_env, ewald_env, atomic_kind_set,& 00241 particle_set, cell, e_rspace, my_task,& 00242 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles,& 00243 forces_glob, pv_glob, efield0_sr, efield1_sr, efield2_sr, error) 00244 energy_glob = energy_glob + e_rspace 00245 00246 IF (do_correction_bonded) THEN 00247 ! The corrections for bonded interactions are stored in the Real Space 00248 ! (Short-Range) part of the fields array. 00249 CALL ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, & 00250 cell, e_bonded, my_task, do_forces, do_efield, do_stress, & 00251 charges, dipoles, quadrupoles, forces_glob, pv_glob, & 00252 efield0_sr, efield1_sr, efield2_sr, error) 00253 energy_glob = energy_glob + e_bonded 00254 END IF 00255 END IF 00256 00257 e_neut = 0.0_dp 00258 e_self = 0.0_dp 00259 energy_local = 0.0_dp 00260 IF (.NOT.debug_r_space) THEN 00261 ! Compute the Reciprocal Space (Long-Range) part of the Ewald sum 00262 CALL ewald_multipole_LR(ewald_env, ewald_pw, cell, particle_set, & 00263 local_particles, energy_local, my_task, do_forces, do_efield, do_stress,& 00264 charges, dipoles, quadrupoles, forces_local, pv_local, efield0_lr, efield1_lr,& 00265 efield2_lr, error) 00266 00267 ! Self-Interactions corrections 00268 CALL ewald_multipole_self (ewald_env, cell, local_particles, e_self, & 00269 e_neut, my_task, do_efield, radii, charges, dipoles, quadrupoles, & 00270 efield0_lr, efield1_lr, efield2_lr, error) 00271 END IF 00272 00273 ! Sumup energy contributions for possible IO 00274 CALL ewald_env_get (ewald_env, group=group, error=error) 00275 energy_glob_t = energy_glob 00276 e_rspace_t = e_rspace 00277 e_bonded_t = e_bonded 00278 CALL mp_sum(energy_glob_t, group) 00279 CALL mp_sum(e_rspace_t, group) 00280 CALL mp_sum(e_bonded_t, group) 00281 ! Print some info about energetics 00282 CALL ewald_multipole_print (iw, energy_local, e_rspace_t, e_bonded_t, e_self, e_neut) 00283 00284 ! Gather the components of the potential, fields and electrostatic field gradients 00285 IF (do_efield) THEN 00286 IF (PRESENT(efield0)) THEN 00287 efield0 = efield0_sr + efield0_lr 00288 CALL mp_sum(efield0, group) 00289 DEALLOCATE (efield0_sr, stat=stat) 00290 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00291 DEALLOCATE (efield0_lr, stat=stat) 00292 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00293 END IF 00294 IF (PRESENT(efield1)) THEN 00295 efield1 = efield1_sr + efield1_lr 00296 CALL mp_sum(efield1, group) 00297 DEALLOCATE (efield1_sr, stat=stat) 00298 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00299 DEALLOCATE (efield1_lr, stat=stat) 00300 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00301 END IF 00302 IF (PRESENT(efield2)) THEN 00303 efield2 = efield2_sr + efield2_lr 00304 CALL mp_sum(efield2, group) 00305 DEALLOCATE (efield2_sr, stat=stat) 00306 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00307 DEALLOCATE (efield2_lr, stat=stat) 00308 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00309 END IF 00310 END IF 00311 CALL timestop(handle) 00312 END SUBROUTINE ewald_multipole_evaluate 00313 00314 ! ***************************************************************************** 00319 SUBROUTINE ewald_multipole_SR (nonbond_env, ewald_env, atomic_kind_set,& 00320 particle_set, cell, energy, task,& 00321 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles,& 00322 forces, pv, efield0, efield1, efield2, error) 00323 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00324 TYPE(ewald_environment_type), POINTER :: ewald_env 00325 TYPE(atomic_kind_type), DIMENSION(:), 00326 OPTIONAL, POINTER :: atomic_kind_set 00327 TYPE(particle_type), POINTER :: particle_set(:) 00328 TYPE(cell_type), POINTER :: cell 00329 REAL(KIND=dp), INTENT(INOUT) :: energy 00330 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 00331 LOGICAL, INTENT(IN) :: do_forces, do_efield, 00332 do_stress 00333 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00334 POINTER :: radii, charges 00335 REAL(KIND=dp), DIMENSION(:, :), 00336 OPTIONAL, POINTER :: dipoles 00337 REAL(KIND=dp), DIMENSION(:, :, :), 00338 OPTIONAL, POINTER :: quadrupoles 00339 REAL(KIND=dp), DIMENSION(:, :), 00340 INTENT(INOUT), OPTIONAL :: forces, pv 00341 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 00342 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 00343 TYPE(cp_error_type), INTENT(inout) :: error 00344 00345 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_SR', 00346 routineP = moduleN//':'//routineN 00347 00348 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, i, iend, igrp, ikind, 00349 ilist, ipair, istart, itype_ij, itype_ji, jkind, k, kind_a, kind_b, kk, 00350 nkdamp_ij, nkdamp_ji, nkinds, npairs 00351 INTEGER, DIMENSION(:, :), POINTER :: list 00352 LOGICAL :: do_efield0, do_efield1, 00353 do_efield2, failure, 00354 force_eval 00355 REAL(KIND=dp) :: alpha, beta, ch_i, ch_j, dampa_ij, dampa_ji, dampaexpi, 00356 dampaexpj, dampfac_ij, dampfac_ji, dampfuncdiffi, dampfuncdiffj, 00357 dampfunci, dampfuncj, dampsumfi, dampsumfj, ef0_i, ef0_j, eloc, fac, 00358 fac_ij, factorial, ir, irab2, ptens11, ptens12, ptens13, ptens21, 00359 ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, rab2_max, radius, 00360 rcut, tij, tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, 00361 tmp31, tmp32, tmp33, tmp_ij, tmp_ji, xf 00362 REAL(KIND=dp), DIMENSION(0:5) :: f 00363 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi, damptij_a, 00364 damptji_a, dp_i, dp_j, ef1_i, 00365 ef1_j, fr, rab, tij_a 00366 REAL(KIND=dp), DIMENSION(3, 3) :: damptij_ab, damptji_ab, 00367 ef2_i, ef2_j, qp_i, qp_j, 00368 tij_ab 00369 REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc 00370 REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd 00371 REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde 00372 TYPE(damping_type) :: damping_ij, damping_ji 00373 TYPE(fist_neighbor_type), POINTER :: nonbonded 00374 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 00375 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, 00376 r_last_update_pbc 00377 00378 failure = .FALSE. 00379 CALL timeset ( routineN, handle ) 00380 NULLIFY(nonbonded,r_last_update, r_last_update_pbc) 00381 do_efield0 = do_efield.AND.ASSOCIATED(efield0) 00382 do_efield1 = do_efield.AND.ASSOCIATED(efield1) 00383 do_efield2 = do_efield.AND.ASSOCIATED(efield2) 00384 IF (do_stress) THEN 00385 ptens11 = 0.0_dp ; ptens12 = 0.0_dp ; ptens13 = 0.0_dp 00386 ptens21 = 0.0_dp ; ptens22 = 0.0_dp ; ptens23 = 0.0_dp 00387 ptens31 = 0.0_dp ; ptens32 = 0.0_dp ; ptens33 = 0.0_dp 00388 END IF 00389 ! Get nonbond_env info 00390 CALL fist_nonbond_env_get (nonbond_env, nonbonded=nonbonded, natom_types = nkinds,& 00391 r_last_update=r_last_update,r_last_update_pbc=r_last_update_pbc, error=error) 00392 CALL ewald_env_get (ewald_env, alpha=alpha, rcut=rcut, error=error) 00393 rab2_max = rcut**2 00394 IF (debug_r_space) THEN 00395 rab2_max = HUGE(0.0_dp) 00396 END IF 00397 ! Starting the force loop 00398 Lists: DO ilist=1,nonbonded%nlists 00399 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00400 npairs=neighbor_kind_pair%npairs 00401 IF (npairs ==0) CYCLE 00402 list => neighbor_kind_pair%list 00403 cvi = neighbor_kind_pair%cell_vector 00404 CALL matvec_3x3(cell_v, cell%hmat, cvi) 00405 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00406 istart = neighbor_kind_pair%grp_kind_start(igrp) 00407 iend = neighbor_kind_pair%grp_kind_end(igrp) 00408 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00409 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00410 00411 itype_ij=no_damping 00412 nkdamp_ij=1 00413 dampa_ij=1.0_dp 00414 dampfac_ij=0.0_dp 00415 00416 itype_ji=no_damping 00417 nkdamp_ji=1 00418 dampa_ji=1.0_dp 00419 dampfac_ji=0.0_dp 00420 IF (PRESENT(atomic_kind_set)) THEN 00421 IF (ASSOCIATED(atomic_kind_set(jkind)%damping)) THEN 00422 damping_ij=atomic_kind_set(jkind)%damping%damp(ikind) 00423 itype_ij=damping_ij%itype 00424 nkdamp_ij=damping_ij%order 00425 dampa_ij=damping_ij%bij 00426 dampfac_ij=damping_ij%cij 00427 END IF 00428 00429 IF (ASSOCIATED(atomic_kind_set(ikind)%damping)) THEN 00430 damping_ji=atomic_kind_set(ikind)%damping%damp(jkind) 00431 itype_ji=damping_ji%itype 00432 nkdamp_ji=damping_ji%order 00433 dampa_ji=damping_ji%bij 00434 dampfac_ji=damping_ji%cij 00435 END IF 00436 END IF 00437 00438 Pairs: DO ipair = istart, iend 00439 IF (ipair <= neighbor_kind_pair%nscale) THEN 00440 ! scale the electrostatic interaction if needed 00441 ! (most often scaled to zero) 00442 fac_ij = neighbor_kind_pair%ei_scale(ipair) 00443 IF (fac_ij<=0) CYCLE 00444 ELSE 00445 fac_ij = 1.0_dp 00446 END IF 00447 atom_a = list(1,ipair) 00448 atom_b = list(2,ipair) 00449 kind_a=particle_set(atom_a)%atomic_kind%kind_number 00450 kind_b=particle_set(atom_b)%atomic_kind%kind_number 00451 IF (atom_a==atom_b) fac_ij = 0.5_dp 00452 rab = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 00453 rab = rab + cell_v 00454 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2 00455 IF (rab2 <= rab2_max) THEN 00456 IF (PRESENT(radii)) THEN 00457 radius = SQRT(radii(atom_a)*radii(atom_a) + radii(atom_b)*radii(atom_b)) 00458 ELSE 00459 radius = 0.0_dp 00460 END IF 00461 IF (radius > 0.0_dp) THEN 00462 beta = sqrthalf/radius 00463 #define SCREENED_COULOMB_GAUSS 00464 #define STORE_ENERGY 00465 #define STORE_FORCES 00466 #include "ewalds_multipole_sr.f90" 00467 #undef STORE_FORCES 00468 #undef STORE_ENERGY 00469 #undef SCREENED_COULOMB_GAUSS 00470 ELSE 00471 #define SCREENED_COULOMB_ERFC 00472 #define STORE_ENERGY 00473 #define STORE_FORCES 00474 #define DAMPING 00475 #include "ewalds_multipole_sr.f90" 00476 #undef DAMPING 00477 #undef STORE_FORCES 00478 #undef STORE_ENERGY 00479 #undef SCREENED_COULOMB_ERFC 00480 END IF 00481 END IF 00482 END DO Pairs 00483 END DO Kind_Group_Loop 00484 END DO Lists 00485 IF (do_stress) THEN 00486 pv(1,1) = pv(1,1) + ptens11 00487 pv(1,2) = pv(1,2) + (ptens12+ptens21)*0.5_dp 00488 pv(1,3) = pv(1,3) + (ptens13+ptens31)*0.5_dp 00489 pv(2,1) = pv(1,2) 00490 pv(2,2) = pv(2,2) + ptens22 00491 pv(2,3) = pv(2,3) + (ptens23+ptens32)*0.5_dp 00492 pv(3,1) = pv(1,3) 00493 pv(3,2) = pv(2,3) 00494 pv(3,3) = pv(3,3) + ptens33 00495 END IF 00496 00497 CALL timestop ( handle ) 00498 END SUBROUTINE ewald_multipole_SR 00499 00500 ! ***************************************************************************** 00505 SUBROUTINE ewald_multipole_bonded (nonbond_env, particle_set, ewald_env, & 00506 cell, energy, task, do_forces, do_efield, do_stress, charges, & 00507 dipoles, quadrupoles, forces, pv, efield0, efield1, efield2, error) 00508 00509 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00510 TYPE(particle_type), POINTER :: particle_set( : ) 00511 TYPE(ewald_environment_type), POINTER :: ewald_env 00512 TYPE(cell_type), POINTER :: cell 00513 REAL(KIND=dp), INTENT(INOUT) :: energy 00514 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 00515 LOGICAL, INTENT(IN) :: do_forces, do_efield, 00516 do_stress 00517 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00518 POINTER :: charges 00519 REAL(KIND=dp), DIMENSION(:, :), 00520 OPTIONAL, POINTER :: dipoles 00521 REAL(KIND=dp), DIMENSION(:, :, :), 00522 OPTIONAL, POINTER :: quadrupoles 00523 REAL(KIND=dp), DIMENSION(:, :), 00524 INTENT(INOUT), OPTIONAL :: forces, pv 00525 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 00526 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 00527 TYPE(cp_error_type), INTENT(inout) :: error 00528 00529 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_bonded', 00530 routineP = moduleN//':'//routineN 00531 00532 INTEGER :: a, atom_a, atom_b, b, c, d, 00533 e, handle, i, iend, igrp, 00534 ilist, ipair, istart, k, 00535 nscale 00536 INTEGER, DIMENSION(:, :), POINTER :: list 00537 LOGICAL :: do_efield0, do_efield1, 00538 do_efield2, force_eval 00539 REAL(KIND=dp) :: alpha, ch_i, ch_j, ef0_i, ef0_j, eloc, fac, fac_ij, ir, 00540 irab2, ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, 00541 ptens32, ptens33, r, rab2, tij, tmp, tmp1, tmp11, tmp12, tmp13, tmp2, 00542 tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, tmp_ji 00543 REAL(KIND=dp), DIMENSION(0:5) :: f 00544 REAL(KIND=dp), DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, 00545 rab, tij_a 00546 REAL(KIND=dp), DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, 00547 tij_ab 00548 REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc 00549 REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd 00550 REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde 00551 TYPE(fist_neighbor_type), POINTER :: nonbonded 00552 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 00553 00554 CALL timeset ( routineN, handle ) 00555 do_efield0 = do_efield.AND.ASSOCIATED(efield0) 00556 do_efield1 = do_efield.AND.ASSOCIATED(efield1) 00557 do_efield2 = do_efield.AND.ASSOCIATED(efield2) 00558 IF (do_stress) THEN 00559 ptens11 = 0.0_dp ; ptens12 = 0.0_dp ; ptens13 = 0.0_dp 00560 ptens21 = 0.0_dp ; ptens22 = 0.0_dp ; ptens23 = 0.0_dp 00561 ptens31 = 0.0_dp ; ptens32 = 0.0_dp ; ptens33 = 0.0_dp 00562 END IF 00563 CALL ewald_env_get (ewald_env, alpha=alpha, error=error) 00564 CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded, error=error) 00565 00566 ! Starting the force loop 00567 Lists: DO ilist=1,nonbonded%nlists 00568 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00569 nscale = neighbor_kind_pair%nscale 00570 IF (nscale==0) CYCLE 00571 list => neighbor_kind_pair%list 00572 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00573 istart = neighbor_kind_pair%grp_kind_start(igrp) 00574 IF (istart>nscale) CYCLE 00575 iend = MIN(neighbor_kind_pair%grp_kind_end(igrp), nscale) 00576 Pairs: DO ipair = istart, iend 00577 ! only use pairs that are (partially) excluded for electrostatics 00578 fac_ij = -1.0_dp + neighbor_kind_pair%ei_scale(ipair) 00579 IF (fac_ij>=0) CYCLE 00580 00581 atom_a = list(1,ipair) 00582 atom_b = list(2,ipair) 00583 00584 rab = particle_set(atom_b)%r - particle_set(atom_a)%r 00585 rab = pbc(rab, cell) 00586 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2 00587 #define SCREENED_COULOMB_ERF 00588 #define STORE_ENERGY 00589 #define STORE_FORCES 00590 #include "ewalds_multipole_sr.f90" 00591 #undef STORE_FORCES 00592 #undef STORE_ENERGY 00593 #undef SCREENED_COULOMB_ERF 00594 END DO Pairs 00595 END DO Kind_Group_Loop 00596 END DO Lists 00597 IF (do_stress) THEN 00598 pv(1,1) = pv(1,1) + ptens11 00599 pv(1,2) = pv(1,2) + (ptens12+ptens21)*0.5_dp 00600 pv(1,3) = pv(1,3) + (ptens13+ptens31)*0.5_dp 00601 pv(2,1) = pv(1,2) 00602 pv(2,2) = pv(2,2) + ptens22 00603 pv(2,3) = pv(2,3) + (ptens23+ptens32)*0.5_dp 00604 pv(3,1) = pv(1,3) 00605 pv(3,2) = pv(2,3) 00606 pv(3,3) = pv(3,3) + ptens33 00607 END IF 00608 00609 CALL timestop ( handle ) 00610 END SUBROUTINE ewald_multipole_bonded 00611 00612 ! ***************************************************************************** 00617 SUBROUTINE ewald_multipole_LR(ewald_env, ewald_pw, cell, particle_set, & 00618 local_particles, energy, task, do_forces, do_efield, do_stress, & 00619 charges, dipoles, quadrupoles, forces, pv, efield0, efield1, efield2,& 00620 error) 00621 TYPE(ewald_environment_type), POINTER :: ewald_env 00622 TYPE(ewald_pw_type), POINTER :: ewald_pw 00623 TYPE(cell_type), POINTER :: cell 00624 TYPE(particle_type), POINTER :: particle_set( : ) 00625 TYPE(distribution_1d_type), POINTER :: local_particles 00626 REAL(KIND=dp), INTENT(INOUT) :: energy 00627 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 00628 LOGICAL, INTENT(IN) :: do_forces, do_efield, 00629 do_stress 00630 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00631 POINTER :: charges 00632 REAL(KIND=dp), DIMENSION(:, :), 00633 OPTIONAL, POINTER :: dipoles 00634 REAL(KIND=dp), DIMENSION(:, :, :), 00635 OPTIONAL, POINTER :: quadrupoles 00636 REAL(KIND=dp), DIMENSION(:, :), 00637 INTENT(INOUT), OPTIONAL :: forces, pv 00638 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 00639 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 00640 TYPE(cp_error_type), INTENT(inout) :: error 00641 00642 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_LR', 00643 routineP = moduleN//':'//routineN 00644 00645 COMPLEX(KIND=dp) :: atm_factor, atm_factor_st(3), 00646 cnjg_fac, fac, summe_tmp 00647 COMPLEX(KIND=dp), ALLOCATABLE, 00648 DIMENSION(:) :: summe_ef 00649 COMPLEX(KIND=dp), ALLOCATABLE, 00650 DIMENSION(:, :) :: summe_st 00651 INTEGER :: gpt, group, handle, iparticle, iparticle_kind, 00652 iparticle_local, lp, mp, nnodes, node, np, nparticle_kind, 00653 nparticle_local, stat 00654 INTEGER, DIMENSION(:, :), POINTER :: bds 00655 LOGICAL :: do_efield0, do_efield1, 00656 do_efield2, failure 00657 REAL(KIND=dp) :: alpha, denom, dipole_t(3), 00658 f0, factor, four_alpha_sq, 00659 gauss, pref, q_t, tmp, trq_t 00660 REAL(KIND=dp), DIMENSION(3) :: tmp_v, vec 00661 REAL(KIND=dp), DIMENSION(3, 3) :: pv_tmp 00662 REAL(KIND=dp), DIMENSION(:, :, :), 00663 POINTER :: rho0 00664 TYPE(dg_rho0_type), POINTER :: dg_rho0 00665 TYPE(dg_type), POINTER :: dg 00666 TYPE(pw_grid_type), POINTER :: pw_grid 00667 TYPE(pw_pool_type), POINTER :: pw_pool 00668 TYPE(structure_factor_type) :: exp_igr 00669 00670 CALL timeset(routineN,handle) 00671 do_efield0 = do_efield.AND.ASSOCIATED(efield0) 00672 do_efield1 = do_efield.AND.ASSOCIATED(efield1) 00673 do_efield2 = do_efield.AND.ASSOCIATED(efield2) 00674 00675 ! Gathering data from the ewald environment 00676 CALL ewald_env_get (ewald_env, alpha=alpha, group=group, error=error) 00677 CALL ewald_pw_get (ewald_pw, pw_big_pool=pw_pool, dg=dg) 00678 CALL dg_get (dg, dg_rho0=dg_rho0) 00679 rho0 => dg_rho0%density%pw%cr3d 00680 pw_grid => pw_pool%pw_grid 00681 bds => pw_grid%bounds 00682 00683 ! Allocation of working arrays 00684 nparticle_kind = SIZE(local_particles%n_el) 00685 nnodes = 0 00686 DO iparticle_kind = 1, nparticle_kind 00687 nnodes = nnodes + local_particles%n_el(iparticle_kind) 00688 ENDDO 00689 CALL structure_factor_allocate(pw_grid%bounds, nnodes, exp_igr) 00690 00691 ALLOCATE (summe_ef(1:pw_grid%ngpts_cut), stat=stat) 00692 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00693 summe_ef = CMPLX (0.0_dp, 0.0_dp,KIND=dp) 00694 ! Stress Tensor 00695 IF (do_stress) THEN 00696 pv_tmp = 0.0_dp 00697 ALLOCATE (summe_st(3,1:pw_grid%ngpts_cut), stat=stat) 00698 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00699 summe_st = CMPLX (0.0_dp, 0.0_dp,KIND=dp) 00700 END IF 00701 00702 ! Defining four_alpha_sq 00703 four_alpha_sq = 4.0_dp * alpha **2 00704 dipole_t = 0.0_dp 00705 q_t = 0.0_dp 00706 trq_t = 0.0_dp 00707 ! Zero node count 00708 node = 0 00709 DO iparticle_kind = 1, nparticle_kind 00710 nparticle_local = local_particles%n_el(iparticle_kind) 00711 DO iparticle_local=1,nparticle_local 00712 node = node + 1 00713 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00714 CALL matvec_3x3 (vec, cell%h_inv, particle_set(iparticle)%r) 00715 CALL structure_factor_evaluate (vec, pw_grid%npts, exp_igr%lb, & 00716 exp_igr%ex(:,node), exp_igr%ey(:,node), exp_igr%ez(:,node)) 00717 00718 ! Computing the total charge, dipole and quadrupole trace (if any) 00719 IF (ANY(task(1,:))) THEN 00720 q_t = q_t + charges( iparticle) 00721 END IF 00722 IF (ANY(task(2,:))) THEN 00723 dipole_t = dipole_t + dipoles(:,iparticle) 00724 END IF 00725 IF (ANY(task(3,:))) THEN 00726 trq_t = trq_t + quadrupoles(1,1,iparticle)+& 00727 quadrupoles(2,2,iparticle)+& 00728 quadrupoles(3,3,iparticle) 00729 END IF 00730 END DO 00731 END DO 00732 00733 ! Looping over the positive g-vectors 00734 DO gpt = 1, pw_grid%ngpts_cut 00735 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt)) 00736 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt)) 00737 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt)) 00738 00739 lp = lp + bds(1,1) 00740 mp = mp + bds(1,2) 00741 np = np + bds(1,3) 00742 00743 ! Initializing sum to be used in the energy and force 00744 node = 0 00745 DO iparticle_kind = 1, nparticle_kind 00746 nparticle_local = local_particles%n_el(iparticle_kind) 00747 DO iparticle_local=1,nparticle_local 00748 node = node + 1 00749 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00750 ! Density for energy and forces 00751 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges,& 00752 dipoles, quadrupoles, error) 00753 summe_tmp = exp_igr%ex(lp,node)*exp_igr%ey(mp,node)*exp_igr%ez(np,node) 00754 summe_ef(gpt) = summe_ef(gpt) + atm_factor*summe_tmp 00755 00756 ! Precompute pseudo-density for stress tensor calculation 00757 IF (do_stress) THEN 00758 CALL get_atom_factor_stress(atm_factor_st, pw_grid, gpt, iparticle, task,& 00759 dipoles, quadrupoles, error) 00760 summe_st(1:3,gpt) = summe_st(1:3,gpt) + atm_factor_st(1:3) *summe_tmp 00761 END IF 00762 END DO 00763 END DO 00764 END DO 00765 CALL mp_sum ( q_t, group) 00766 CALL mp_sum (dipole_t, group) 00767 CALL mp_sum ( trq_t, group) 00768 CALL mp_sum (summe_ef, group) 00769 IF (do_stress) CALL mp_sum (summe_st, group) 00770 00771 ! Looping over the positive g-vectors 00772 DO gpt = 1, pw_grid%ngpts_cut 00773 ! computing the potential energy 00774 lp = pw_grid%mapl%pos(pw_grid%g_hat(1,gpt)) 00775 mp = pw_grid%mapm%pos(pw_grid%g_hat(2,gpt)) 00776 np = pw_grid%mapn%pos(pw_grid%g_hat(3,gpt)) 00777 00778 lp = lp + bds(1,1) 00779 mp = mp + bds(1,2) 00780 np = np + bds(1,3) 00781 00782 IF (pw_grid%gsq(gpt) == 0.0_dp) THEN 00783 ! G=0 vector for dipole-dipole and charge-quadrupole 00784 energy = energy + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t,dipole_t)& 00785 - (1.0_dp/9.0_dp)*q_t*trq_t 00786 ! Stress tensor 00787 IF (do_stress) THEN 00788 pv_tmp(1,1) = pv_tmp(1,1) + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t,dipole_t) 00789 pv_tmp(2,2) = pv_tmp(2,2) + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t,dipole_t) 00790 pv_tmp(3,3) = pv_tmp(3,3) + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t,dipole_t) 00791 END IF 00792 ! Corrections for G=0 to potential, field and field gradient 00793 IF (do_efield.AND.(debug_e_field_en.OR.(.NOT.debug_this_module))) THEN 00794 ! This term is important and may give problems if one is debugging 00795 ! VS finite differences since it comes from a residual integral in 00796 ! the complex plane (cannot be reproduced with finite differences) 00797 node = 0 00798 DO iparticle_kind = 1, nparticle_kind 00799 nparticle_local = local_particles%n_el(iparticle_kind) 00800 DO iparticle_local=1,nparticle_local 00801 node = node + 1 00802 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00803 00804 ! Potential 00805 IF (do_efield0) THEN 00806 efield0( iparticle) = efield0( iparticle) 00807 END IF 00808 ! Electrostatic field 00809 IF (do_efield1) THEN 00810 efield1(1:3,iparticle) = efield1(1:3,iparticle) - (1.0_dp/6.0_dp)*dipole_t(1:3) 00811 END IF 00812 ! Electrostatic field gradients 00813 IF (do_efield2) THEN 00814 efield2(1,iparticle) = efield2(1,iparticle) - (1.0_dp/(18.0_dp))*q_t 00815 efield2(5,iparticle) = efield2(5,iparticle) - (1.0_dp/(18.0_dp))*q_t 00816 efield2(9,iparticle) = efield2(9,iparticle) - (1.0_dp/(18.0_dp))*q_t 00817 END IF 00818 END DO 00819 END DO 00820 END IF 00821 CYCLE 00822 END IF 00823 gauss = (rho0(lp,mp,np) * pw_grid%vol)**2 / pw_grid%gsq(gpt) 00824 factor = gauss * REAL(summe_ef(gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00825 energy = energy + factor 00826 00827 IF (do_forces.OR.do_efield) THEN 00828 node = 0 00829 DO iparticle_kind = 1, nparticle_kind 00830 nparticle_local = local_particles%n_el(iparticle_kind) 00831 DO iparticle_local=1,nparticle_local 00832 node = node + 1 00833 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00834 fac = exp_igr%ex(lp,node)*exp_igr%ey(mp,node)*exp_igr%ez(np,node) 00835 cnjg_fac = CONJG(fac) 00836 00837 ! Forces 00838 IF (do_forces) THEN 00839 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges,& 00840 dipoles, quadrupoles, error) 00841 00842 tmp = gauss * AIMAG(summe_ef(gpt) * (cnjg_fac * CONJG(atm_factor))) 00843 forces(1,node) = forces(1,node) + tmp * pw_grid%g(1,gpt) 00844 forces(2,node) = forces(2,node) + tmp * pw_grid%g(2,gpt) 00845 forces(3,node) = forces(3,node) + tmp * pw_grid%g(3,gpt) 00846 END IF 00847 00848 ! Electric field 00849 IF (do_efield) THEN 00850 ! Potential 00851 IF (do_efield0) THEN 00852 efield0( iparticle) = efield0( iparticle) + gauss * REAL(fac*CONJG(summe_ef(gpt)),KIND=dp) 00853 END IF 00854 ! Electric field 00855 IF (do_efield1) THEN 00856 tmp = AIMAG(fac*CONJG(summe_ef(gpt)))*gauss 00857 efield1(1,iparticle) = efield1(1,iparticle) - tmp * pw_grid%g(1,gpt) 00858 efield1(2,iparticle) = efield1(2,iparticle) - tmp * pw_grid%g(2,gpt) 00859 efield1(3,iparticle) = efield1(3,iparticle) - tmp * pw_grid%g(3,gpt) 00860 END IF 00861 ! Electric field gradient 00862 IF (do_efield2) THEN 00863 tmp_v(1) = REAL(fac*CONJG(summe_ef(gpt)),KIND=dp)*pw_grid%g(1,gpt)*gauss 00864 tmp_v(2) = REAL(fac*CONJG(summe_ef(gpt)),KIND=dp)*pw_grid%g(2,gpt)*gauss 00865 tmp_v(3) = REAL(fac*CONJG(summe_ef(gpt)),KIND=dp)*pw_grid%g(3,gpt)*gauss 00866 00867 efield2(1,iparticle) = efield2(1,iparticle) + tmp_v(1) * pw_grid%g(1,gpt) 00868 efield2(2,iparticle) = efield2(2,iparticle) + tmp_v(1) * pw_grid%g(2,gpt) 00869 efield2(3,iparticle) = efield2(3,iparticle) + tmp_v(1) * pw_grid%g(3,gpt) 00870 efield2(4,iparticle) = efield2(4,iparticle) + tmp_v(2) * pw_grid%g(1,gpt) 00871 efield2(5,iparticle) = efield2(5,iparticle) + tmp_v(2) * pw_grid%g(2,gpt) 00872 efield2(6,iparticle) = efield2(6,iparticle) + tmp_v(2) * pw_grid%g(3,gpt) 00873 efield2(7,iparticle) = efield2(7,iparticle) + tmp_v(3) * pw_grid%g(1,gpt) 00874 efield2(8,iparticle) = efield2(8,iparticle) + tmp_v(3) * pw_grid%g(2,gpt) 00875 efield2(9,iparticle) = efield2(9,iparticle) + tmp_v(3) * pw_grid%g(3,gpt) 00876 END IF 00877 END IF 00878 END DO 00879 END DO 00880 END IF 00881 00882 ! Compute the virial P*V 00883 IF (do_stress) THEN 00884 ! The Stress Tensor can be decomposed in two main components. 00885 ! The first one is just a normal ewald component for reciprocal space 00886 denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt) 00887 pv_tmp(1,1) = pv_tmp(1,1) + factor*(1.0_dp - 2.0_dp*pw_grid%g(1,gpt)*pw_grid%g(1,gpt)*denom) 00888 pv_tmp(1,2) = pv_tmp(1,2) - factor*( 2.0_dp*pw_grid%g(1,gpt)*pw_grid%g(2,gpt)*denom) 00889 pv_tmp(1,3) = pv_tmp(1,3) - factor*( 2.0_dp*pw_grid%g(1,gpt)*pw_grid%g(3,gpt)*denom) 00890 pv_tmp(2,1) = pv_tmp(2,1) - factor*( 2.0_dp*pw_grid%g(2,gpt)*pw_grid%g(1,gpt)*denom) 00891 pv_tmp(2,2) = pv_tmp(2,2) + factor*(1.0_dp - 2.0_dp*pw_grid%g(2,gpt)*pw_grid%g(2,gpt)*denom) 00892 pv_tmp(2,3) = pv_tmp(2,3) - factor*( 2.0_dp*pw_grid%g(2,gpt)*pw_grid%g(3,gpt)*denom) 00893 pv_tmp(3,1) = pv_tmp(3,1) - factor*( 2.0_dp*pw_grid%g(3,gpt)*pw_grid%g(1,gpt)*denom) 00894 pv_tmp(3,2) = pv_tmp(3,2) - factor*( 2.0_dp*pw_grid%g(3,gpt)*pw_grid%g(2,gpt)*denom) 00895 pv_tmp(3,3) = pv_tmp(3,3) + factor*(1.0_dp - 2.0_dp*pw_grid%g(3,gpt)*pw_grid%g(3,gpt)*denom) 00896 ! The second one can be written in the following way 00897 f0 = 2.0_dp * gauss 00898 pv_tmp(1,1) = pv_tmp(1,1) + f0 * pw_grid%g(1,gpt) * REAL(summe_st(1,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00899 pv_tmp(1,2) = pv_tmp(1,2) + f0 * pw_grid%g(1,gpt) * REAL(summe_st(2,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00900 pv_tmp(1,3) = pv_tmp(1,3) + f0 * pw_grid%g(1,gpt) * REAL(summe_st(3,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00901 pv_tmp(2,1) = pv_tmp(2,1) + f0 * pw_grid%g(2,gpt) * REAL(summe_st(1,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00902 pv_tmp(2,2) = pv_tmp(2,2) + f0 * pw_grid%g(2,gpt) * REAL(summe_st(2,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00903 pv_tmp(2,3) = pv_tmp(2,3) + f0 * pw_grid%g(2,gpt) * REAL(summe_st(3,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00904 pv_tmp(3,1) = pv_tmp(3,1) + f0 * pw_grid%g(3,gpt) * REAL(summe_st(1,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00905 pv_tmp(3,2) = pv_tmp(3,2) + f0 * pw_grid%g(3,gpt) * REAL(summe_st(2,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00906 pv_tmp(3,3) = pv_tmp(3,3) + f0 * pw_grid%g(3,gpt) * REAL(summe_st(3,gpt) * CONJG(summe_ef(gpt)),KIND=dp) 00907 END IF 00908 END DO 00909 pref = fourpi/pw_grid%vol 00910 energy = energy * pref 00911 00912 CALL structure_factor_deallocate (exp_igr) 00913 DEALLOCATE (summe_ef, stat=stat) 00914 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00915 IF (do_stress) THEN 00916 pv_tmp = pv_tmp * pref 00917 ! Symmetrize the tensor 00918 pv(1,1) = pv(1,1) + pv_tmp(1,1) 00919 pv(1,2) = pv(1,2) + (pv_tmp(1,2) + pv_tmp(2,1))*0.5_dp 00920 pv(1,3) = pv(1,3) + (pv_tmp(1,3) + pv_tmp(3,1))*0.5_dp 00921 pv(2,1) = pv(1,2) 00922 pv(2,2) = pv(2,2) + pv_tmp(2,2) 00923 pv(2,3) = pv(2,3) + (pv_tmp(2,3) + pv_tmp(3,2))*0.5_dp 00924 pv(3,1) = pv(1,3) 00925 pv(3,2) = pv(2,3) 00926 pv(3,3) = pv(3,3) + pv_tmp(3,3) 00927 DEALLOCATE (summe_st, stat=stat) 00928 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00929 END IF 00930 IF (do_forces) THEN 00931 forces = 2.0_dp * forces * pref 00932 END IF 00933 IF (do_efield0) THEN 00934 efield0 = 2.0_dp * efield0 * pref 00935 END IF 00936 IF (do_efield1) THEN 00937 efield1 = 2.0_dp * efield1 * pref 00938 END IF 00939 IF (do_efield2) THEN 00940 efield2 = 2.0_dp * efield2 * pref 00941 END IF 00942 CALL timestop(handle) 00943 00944 END SUBROUTINE ewald_multipole_LR 00945 00946 ! ***************************************************************************** 00952 SUBROUTINE get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges,& 00953 dipoles, quadrupoles, error) 00954 COMPLEX(KIND=dp), INTENT(OUT) :: atm_factor 00955 TYPE(pw_grid_type), POINTER :: pw_grid 00956 INTEGER, INTENT(IN) :: gpt 00957 INTEGER :: iparticle 00958 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 00959 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00960 POINTER :: charges 00961 REAL(KIND=dp), DIMENSION(:, :), 00962 OPTIONAL, POINTER :: dipoles 00963 REAL(KIND=dp), DIMENSION(:, :, :), 00964 OPTIONAL, POINTER :: quadrupoles 00965 TYPE(cp_error_type), INTENT(inout) :: error 00966 00967 CHARACTER(len=*), PARAMETER :: routineN = 'get_atom_factor', 00968 routineP = moduleN//':'//routineN 00969 00970 COMPLEX(KIND=dp) :: tmp 00971 INTEGER :: i, j 00972 00973 atm_factor = CMPLX (0.0_dp, 0.0_dp,KIND=dp) 00974 IF (task(1,1)) THEN 00975 ! Charge 00976 atm_factor = atm_factor + charges(iparticle) 00977 END IF 00978 IF (task(2,2)) THEN 00979 ! Dipole 00980 tmp = CMPLX (0.0_dp, 0.0_dp,KIND=dp) 00981 DO i = 1,3 00982 tmp = tmp + dipoles(i,iparticle)*pw_grid%g(i,gpt) 00983 END DO 00984 atm_factor = atm_factor + tmp * CMPLX(0.0_dp, -1.0_dp, KIND=dp) 00985 END IF 00986 IF (task(3,3)) THEN 00987 ! Quadrupole 00988 tmp = CMPLX (0.0_dp, 0.0_dp,KIND=dp) 00989 DO i = 1,3 00990 DO j = 1,3 00991 tmp = tmp + quadrupoles(j,i,iparticle)*pw_grid%g(j,gpt)*pw_grid%g(i,gpt) 00992 END DO 00993 END DO 00994 atm_factor = atm_factor - 1.0_dp/3.0_dp * tmp 00995 END IF 00996 00997 END SUBROUTINE get_atom_factor 00998 00999 ! ***************************************************************************** 01005 SUBROUTINE get_atom_factor_stress(atm_factor, pw_grid, gpt, iparticle, task,& 01006 dipoles, quadrupoles, error) 01007 COMPLEX(KIND=dp), INTENT(OUT) :: atm_factor(3) 01008 TYPE(pw_grid_type), POINTER :: pw_grid 01009 INTEGER, INTENT(IN) :: gpt 01010 INTEGER :: iparticle 01011 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 01012 REAL(KIND=dp), DIMENSION(:, :), 01013 OPTIONAL, POINTER :: dipoles 01014 REAL(KIND=dp), DIMENSION(:, :, :), 01015 OPTIONAL, POINTER :: quadrupoles 01016 TYPE(cp_error_type), INTENT(inout) :: error 01017 01018 CHARACTER(len=*), PARAMETER :: routineN = 'get_atom_factor_stress', 01019 routineP = moduleN//':'//routineN 01020 01021 INTEGER :: i 01022 01023 atm_factor = CMPLX (0.0_dp, 0.0_dp,KIND=dp) 01024 IF (ANY(task(2,:))) THEN 01025 ! Dipole 01026 atm_factor = dipoles(:,iparticle) * CMPLX(0.0_dp, -1.0_dp, KIND=dp) 01027 END IF 01028 IF (ANY(task(3,:))) THEN 01029 ! Quadrupole 01030 DO i = 1,3 01031 atm_factor(1) = atm_factor(1) - 1.0_dp/3.0_dp *& 01032 (quadrupoles(1,i,iparticle)*pw_grid%g(i,gpt)+& 01033 quadrupoles(i,1,iparticle)*pw_grid%g(i,gpt)) 01034 atm_factor(2) = atm_factor(2) - 1.0_dp/3.0_dp *& 01035 (quadrupoles(2,i,iparticle)*pw_grid%g(i,gpt)+& 01036 quadrupoles(i,2,iparticle)*pw_grid%g(i,gpt)) 01037 atm_factor(3) = atm_factor(3) - 1.0_dp/3.0_dp *& 01038 (quadrupoles(3,i,iparticle)*pw_grid%g(i,gpt)+& 01039 quadrupoles(i,3,iparticle)*pw_grid%g(i,gpt)) 01040 END DO 01041 END IF 01042 01043 END SUBROUTINE get_atom_factor_stress 01044 01045 ! ***************************************************************************** 01050 SUBROUTINE ewald_multipole_self (ewald_env, cell, local_particles, e_self, & 01051 e_neut, task, do_efield, radii, charges, dipoles, quadrupoles, efield0, & 01052 efield1, efield2, error) 01053 TYPE(ewald_environment_type), POINTER :: ewald_env 01054 TYPE(cell_type), INTENT(IN) :: cell 01055 TYPE(distribution_1d_type), POINTER :: local_particles 01056 REAL(KIND=dp), INTENT(OUT) :: e_self, e_neut 01057 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 01058 LOGICAL, INTENT(IN) :: do_efield 01059 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 01060 POINTER :: radii, charges 01061 REAL(KIND=dp), DIMENSION(:, :), 01062 OPTIONAL, POINTER :: dipoles 01063 REAL(KIND=dp), DIMENSION(:, :, :), 01064 OPTIONAL, POINTER :: quadrupoles 01065 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 01066 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 01067 TYPE(cp_error_type), INTENT(inout) :: error 01068 01069 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_self', 01070 routineP = moduleN//':'//routineN 01071 REAL(KIND=dp), PARAMETER :: f23 = 2.0_dp/3.0_dp, 01072 f415 = 4.0_dp/15.0_dp 01073 01074 INTEGER :: ewald_type, group, i, 01075 iparticle, iparticle_kind, 01076 iparticle_local, j, 01077 nparticle_local 01078 LOGICAL :: do_efield0, do_efield1, 01079 do_efield2, lradii 01080 REAL(KIND=dp) :: alpha, ch_qu_self, ch_qu_self_tmp, dipole_self, fac1, 01081 fac2, fac3, fac4, q, q_neutg, q_self, q_sum, qu_qu_self, radius 01082 01083 CALL ewald_env_get ( ewald_env, ewald_type=ewald_type, alpha=alpha,& 01084 group=group, error=error) 01085 01086 do_efield0 = do_efield.AND.ASSOCIATED(efield0) 01087 do_efield1 = do_efield.AND.ASSOCIATED(efield1) 01088 do_efield2 = do_efield.AND.ASSOCIATED(efield2) 01089 q_self = 0.0_dp 01090 q_sum = 0.0_dp 01091 dipole_self = 0.0_dp 01092 ch_qu_self = 0.0_dp 01093 qu_qu_self = 0.0_dp 01094 fac1 = 2.0_dp*alpha*oorootpi 01095 fac2 = 6.0_dp*(f23**2)*(alpha**3)*oorootpi 01096 fac3 = (2.0_dp*oorootpi)*f23*alpha**3 01097 fac4 = (4.0_dp*oorootpi)*f415*alpha**5 01098 lradii = PRESENT(radii) 01099 radius = 0.0_dp 01100 q_neutg = 0.0_dp 01101 DO iparticle_kind=1,SIZE(local_particles%n_el) 01102 nparticle_local = local_particles%n_el(iparticle_kind) 01103 DO iparticle_local=1,nparticle_local 01104 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 01105 IF (ANY(task(1,:))) THEN 01106 ! Charge - Charge 01107 q = charges(iparticle) 01108 IF (lradii) radius = radii(iparticle) 01109 IF (radius>0) THEN 01110 q_neutg = q_neutg + 2.0_dp*q*radius**2 01111 END IF 01112 q_self = q_self + q*q 01113 q_sum = q_sum + q 01114 ! Potential 01115 IF (do_efield0) THEN 01116 efield0( iparticle) = efield0( iparticle) - q*fac1 01117 END IF 01118 01119 IF (task(1,3)) THEN 01120 ! Charge - Quadrupole 01121 ch_qu_self_tmp = 0.0_dp 01122 DO i = 1,3 01123 ch_qu_self_tmp = ch_qu_self_tmp + quadrupoles(i,i,iparticle) * q 01124 END DO 01125 ch_qu_self = ch_qu_self + ch_qu_self_tmp 01126 ! Electric Field Gradient 01127 IF (do_efield2) THEN 01128 efield2(1,iparticle) = efield2(1,iparticle) + fac2*q 01129 efield2(5,iparticle) = efield2(5,iparticle) + fac2*q 01130 efield2(9,iparticle) = efield2(9,iparticle) + fac2*q 01131 END IF 01132 END IF 01133 END IF 01134 IF (ANY(task(2,:))) THEN 01135 ! Dipole - Dipole 01136 DO i = 1,3 01137 dipole_self = dipole_self + dipoles(i,iparticle)**2 01138 END DO 01139 ! Electric Field 01140 IF (do_efield1) THEN 01141 efield1(1,iparticle) = efield1(1,iparticle) + fac3 * dipoles(1,iparticle) 01142 efield1(2,iparticle) = efield1(2,iparticle) + fac3 * dipoles(2,iparticle) 01143 efield1(3,iparticle) = efield1(3,iparticle) + fac3 * dipoles(3,iparticle) 01144 END IF 01145 END IF 01146 IF (ANY(task(3,:))) THEN 01147 ! Quadrupole - Quadrupole 01148 DO i = 1,3 01149 DO j = 1,3 01150 qu_qu_self = qu_qu_self + quadrupoles(j,i,iparticle)**2 01151 END DO 01152 END DO 01153 ! Electric Field Gradient 01154 IF (do_efield2) THEN 01155 efield2(1,iparticle) = efield2(1,iparticle) + fac4 * quadrupoles(1,1,iparticle) 01156 efield2(2,iparticle) = efield2(2,iparticle) + fac4 * quadrupoles(2,1,iparticle) 01157 efield2(3,iparticle) = efield2(3,iparticle) + fac4 * quadrupoles(3,1,iparticle) 01158 efield2(4,iparticle) = efield2(4,iparticle) + fac4 * quadrupoles(1,2,iparticle) 01159 efield2(5,iparticle) = efield2(5,iparticle) + fac4 * quadrupoles(2,2,iparticle) 01160 efield2(6,iparticle) = efield2(6,iparticle) + fac4 * quadrupoles(3,2,iparticle) 01161 efield2(7,iparticle) = efield2(7,iparticle) + fac4 * quadrupoles(1,3,iparticle) 01162 efield2(8,iparticle) = efield2(8,iparticle) + fac4 * quadrupoles(2,3,iparticle) 01163 efield2(9,iparticle) = efield2(9,iparticle) + fac4 * quadrupoles(3,3,iparticle) 01164 END IF 01165 END IF 01166 END DO 01167 END DO 01168 01169 CALL mp_sum (q_neutg, group) 01170 CALL mp_sum (q_self, group) 01171 CALL mp_sum (q_sum, group) 01172 CALL mp_sum (dipole_self, group) 01173 CALL mp_sum (ch_qu_self, group) 01174 CALL mp_sum (qu_qu_self, group) 01175 01176 e_self = -(q_self+f23*(dipole_self-f23*ch_qu_self+f415*qu_qu_self*alpha**2)*alpha**2)*alpha*oorootpi 01177 fac1 = pi/(2.0_dp*cell%deth) 01178 e_neut = -q_sum*fac1*(q_sum/alpha**2 - q_neutg) 01179 01180 ! Correcting Potential for the neutralizing background charge 01181 DO iparticle_kind=1,SIZE(local_particles%n_el) 01182 nparticle_local = local_particles%n_el(iparticle_kind) 01183 DO iparticle_local=1,nparticle_local 01184 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 01185 IF (ANY(task(1,:))) THEN 01186 ! Potential energy 01187 IF (do_efield0) THEN 01188 efield0( iparticle) = efield0( iparticle) - q_sum*2.0_dp*fac1/alpha**2 01189 IF (lradii) radius = radii(iparticle) 01190 IF (radius>0) THEN 01191 q = charges(iparticle) 01192 efield0( iparticle) = efield0( iparticle) + fac1*radius**2*(q_sum+q) 01193 END IF 01194 END IF 01195 END IF 01196 END DO 01197 END DO 01198 END SUBROUTINE ewald_multipole_self 01199 01200 ! ***************************************************************************** 01203 SUBROUTINE ewald_multipole_print(iw, e_gspace, e_rspace, e_bonded, e_self, e_neut) 01204 01205 INTEGER, INTENT(IN) :: iw 01206 REAL(KIND=dp), INTENT(IN) :: e_gspace, e_rspace, e_bonded, 01207 e_self, e_neut 01208 01209 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_print', 01210 routineP = moduleN//':'//routineN 01211 01212 IF (iw>0) THEN 01213 WRITE ( iw, '( A, A )' ) ' *********************************', & 01214 '**********************************************' 01215 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' INITIAL GSPACE ENERGY', & 01216 '[hartree]', '= ', e_gspace 01217 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' INITIAL RSPACE ENERGY', & 01218 '[hartree]', '= ', e_rspace 01219 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' BONDED CORRECTION', & 01220 '[hartree]', '= ', e_bonded 01221 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' SELF ENERGY CORRECTION', & 01222 '[hartree]', '= ', e_self 01223 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' NEUTRALIZ. BCKGR. ENERGY', & 01224 '[hartree]', '= ', e_neut 01225 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' TOTAL ELECTROSTATIC EN.', & 01226 '[hartree]', '= ', e_rspace+e_bonded+e_gspace+e_self+e_neut 01227 WRITE ( iw, '( A, A )' ) ' *********************************', & 01228 '**********************************************' 01229 END IF 01230 END SUBROUTINE ewald_multipole_print 01231 01232 END MODULE ewalds_multipole
1.7.3