|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00011 MODULE ewalds 00012 00013 USE atomic_kind_types, ONLY: atomic_kind_type,& 00014 get_atomic_kind 00015 USE bibliography, ONLY: Ewald1921,& 00016 cite_reference 00017 USE cell_types, ONLY: cell_type 00018 USE dg_rho0_types, ONLY: dg_rho0_type 00019 USE dg_types, ONLY: dg_get,& 00020 dg_type 00021 USE distribution_1d_types, ONLY: distribution_1d_type 00022 USE ewald_environment_types, ONLY: ewald_env_get,& 00023 ewald_environment_type 00024 USE ewald_pw_types, ONLY: ewald_pw_get,& 00025 ewald_pw_type 00026 USE f77_blas 00027 USE input_constants, ONLY: do_ewald_none 00028 USE kinds, ONLY: dp 00029 USE mathconstants, ONLY: fourpi,& 00030 oorootpi,& 00031 pi 00032 USE mathlib, ONLY: matvec_3x3 00033 USE message_passing, ONLY: mp_sum 00034 USE particle_types, ONLY: particle_type 00035 USE pw_grid_types, ONLY: pw_grid_type 00036 USE pw_pool_types, ONLY: pw_pool_type 00037 USE shell_potential_types, ONLY: get_shell,& 00038 shell_kind_type 00039 USE structure_factor_types, ONLY: structure_factor_type 00040 USE structure_factors, ONLY: structure_factor_allocate,& 00041 structure_factor_deallocate,& 00042 structure_factor_evaluate 00043 USE timings, ONLY: timeset,& 00044 timestop 00045 #include "cp_common_uses.h" 00046 00047 IMPLICIT NONE 00048 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00049 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewalds' 00050 00051 PRIVATE 00052 PUBLIC :: ewald_evaluate, ewald_self, ewald_self_atom, ewald_print 00053 00054 CONTAINS 00055 00056 ! ***************************************************************************** 00065 SUBROUTINE ewald_evaluate (ewald_env, ewald_pw, cell, atomic_kind_set, particle_set,& 00066 local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial, charges, e_coulomb, error ) 00067 TYPE(ewald_environment_type), POINTER :: ewald_env 00068 TYPE(ewald_pw_type), POINTER :: ewald_pw 00069 TYPE(cell_type), POINTER :: cell 00070 TYPE(atomic_kind_type), POINTER :: atomic_kind_set( : ) 00071 TYPE(particle_type), POINTER :: particle_set( : ) 00072 TYPE(distribution_1d_type), POINTER :: local_particles 00073 REAL(KIND=dp), DIMENSION(:, :), 00074 INTENT(OUT) :: fg_coulomb 00075 REAL(KIND=dp), INTENT(OUT) :: vg_coulomb 00076 REAL(KIND=dp), DIMENSION(:, :), 00077 INTENT(OUT) :: pv_g 00078 LOGICAL, INTENT(IN) :: use_virial 00079 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00080 POINTER :: charges, e_coulomb 00081 TYPE(cp_error_type), INTENT(inout) :: error 00082 00083 CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_evaluate', 00084 routineP = moduleN//':'//routineN 00085 00086 COMPLEX(KIND=dp) :: snode 00087 COMPLEX(KIND=dp), ALLOCATABLE, 00088 DIMENSION(:) :: summe 00089 INTEGER :: gpt, group, handle, iparticle, iparticle_kind, 00090 iparticle_local, isos, lp, mp, nnodes, node, np, nparticle_kind, 00091 nparticle_local 00092 INTEGER, DIMENSION(:, :), POINTER :: bds 00093 LOGICAL :: atenergy, failure, 00094 use_charge_array 00095 REAL(KIND=dp) :: alpha, denom, e_igdotr, 00096 factor, four_alpha_sq, gauss, 00097 pref, q 00098 REAL(KIND=dp), DIMENSION(3) :: vec 00099 REAL(KIND=dp), DIMENSION(:), POINTER :: charge 00100 REAL(KIND=dp), DIMENSION(:, :, :), 00101 POINTER :: rho0 00102 TYPE(atomic_kind_type), POINTER :: atomic_kind 00103 TYPE(dg_rho0_type), POINTER :: dg_rho0 00104 TYPE(dg_type), POINTER :: dg 00105 TYPE(pw_grid_type), POINTER :: pw_grid 00106 TYPE(pw_pool_type), POINTER :: pw_pool 00107 TYPE(structure_factor_type) :: exp_igr 00108 00109 CALL timeset(routineN,handle) 00110 CALL cite_reference(Ewald1921) 00111 use_charge_array=PRESENT(charges) 00112 IF (use_charge_array) use_charge_array=ASSOCIATED(charges) 00113 atenergy=PRESENT(e_coulomb) 00114 IF (atenergy) atenergy=ASSOCIATED(e_coulomb) 00115 IF (atenergy) e_coulomb=0._dp 00116 00117 ! pointing 00118 CALL ewald_env_get (ewald_env, alpha=alpha, group = group ,error=error) 00119 CALL ewald_pw_get (ewald_pw, pw_big_pool=pw_pool, dg = dg ) 00120 CALL dg_get ( dg, dg_rho0=dg_rho0 ) 00121 rho0 => dg_rho0 % density % pw % cr3d 00122 pw_grid => pw_pool % pw_grid 00123 bds => pw_grid % bounds 00124 00125 ! allocating 00126 nparticle_kind = SIZE ( atomic_kind_set ) 00127 nnodes = 0 00128 DO iparticle_kind = 1, nparticle_kind 00129 nnodes = nnodes + local_particles%n_el(iparticle_kind) 00130 ENDDO 00131 00132 CALL structure_factor_allocate ( pw_grid % bounds, nnodes, exp_igr) 00133 00134 ALLOCATE (summe(1:pw_grid%ngpts_cut),STAT=isos) 00135 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00136 ALLOCATE (charge(1:nnodes),STAT=isos) 00137 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00138 00139 ! Initializing vg_coulomb and fg_coulomb 00140 vg_coulomb = 0.0_dp 00141 fg_coulomb = 0.0_dp 00142 IF (use_virial) pv_g = 0.0_dp 00143 ! defining four_alpha_sq 00144 four_alpha_sq = 4.0_dp * alpha ** 2 00145 ! zero node count 00146 node = 0 00147 DO iparticle_kind = 1, nparticle_kind 00148 nparticle_local = local_particles%n_el(iparticle_kind) 00149 IF (use_charge_array) THEN 00150 DO iparticle_local=1,nparticle_local 00151 node = node + 1 00152 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00153 charge(node) = charges(iparticle) 00154 CALL matvec_3x3 (vec, cell % h_inv, particle_set ( iparticle ) % r ) 00155 CALL structure_factor_evaluate ( vec, pw_grid % npts, exp_igr % lb, & 00156 exp_igr % ex(:,node), exp_igr % ey(:,node), exp_igr % ez(:,node) ) 00157 END DO 00158 ELSE 00159 atomic_kind => atomic_kind_set(iparticle_kind) 00160 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q) 00161 DO iparticle_local=1,nparticle_local 00162 node = node + 1 00163 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 00164 charge(node) = q 00165 CALL matvec_3x3 (vec, cell % h_inv, particle_set ( iparticle ) % r ) 00166 CALL structure_factor_evaluate ( vec, pw_grid % npts, exp_igr % lb, & 00167 exp_igr % ex(:,node), exp_igr % ey(:,node), exp_igr % ez(:,node) ) 00168 END DO 00169 END IF 00170 END DO 00171 00172 summe ( : ) = CMPLX ( 0.0_dp, 0.0_dp,KIND=dp) 00173 ! looping over the positive g-vectors 00174 DO gpt = 1, pw_grid % ngpts_cut 00175 00176 lp = pw_grid % mapl % pos ( pw_grid % g_hat ( 1, gpt ) ) 00177 mp = pw_grid % mapm % pos ( pw_grid % g_hat ( 2, gpt ) ) 00178 np = pw_grid % mapn % pos ( pw_grid % g_hat ( 3, gpt ) ) 00179 00180 lp = lp + bds ( 1, 1 ) 00181 mp = mp + bds ( 1, 2 ) 00182 np = np + bds ( 1, 3 ) 00183 00184 ! initializing sum to be used in the energy and force 00185 DO node = 1, nnodes 00186 summe ( gpt ) = summe ( gpt ) + charge ( node ) * & 00187 ( exp_igr % ex ( lp, node ) & 00188 * exp_igr % ey ( mp, node ) & 00189 * exp_igr % ez ( np, node ) ) 00190 END DO 00191 END DO 00192 CALL mp_sum ( summe, group ) 00193 00194 pref = fourpi / pw_grid % vol 00195 00196 ! looping over the positive g-vectors 00197 DO gpt = 1, pw_grid % ngpts_cut 00198 ! computing the potential energy 00199 lp = pw_grid % mapl % pos ( pw_grid % g_hat ( 1, gpt ) ) 00200 mp = pw_grid % mapm % pos ( pw_grid % g_hat ( 2, gpt ) ) 00201 np = pw_grid % mapn % pos ( pw_grid % g_hat ( 3, gpt ) ) 00202 00203 lp = lp + bds ( 1, 1 ) 00204 mp = mp + bds ( 1, 2 ) 00205 np = np + bds ( 1, 3 ) 00206 00207 IF ( pw_grid % gsq ( gpt ) <= 1.0E-10_dp ) CYCLE 00208 00209 gauss = ( rho0 ( lp, mp, np ) * pw_grid % vol ) ** 2 / pw_grid % gsq ( gpt ) 00210 factor = gauss * REAL ( summe ( gpt ) * CONJG ( summe ( gpt ) ) ) 00211 vg_coulomb = vg_coulomb + factor 00212 00213 ! atomic energies 00214 IF (atenergy) THEN 00215 DO node = 1, nnodes 00216 snode = CONJG ( exp_igr % ex ( lp, node ) & 00217 * exp_igr % ey ( mp, node ) & 00218 * exp_igr % ez ( np, node ) ) 00219 e_coulomb(node) = e_coulomb(node) + gauss*charge(node)*REAL(summe(gpt)*snode) 00220 END DO 00221 END IF 00222 00223 ! computing the force 00224 node = 0 00225 DO node = 1, nnodes 00226 e_igdotr = AIMAG ( summe(gpt) * CONJG & 00227 ( exp_igr % ex ( lp, node ) & 00228 * exp_igr % ey ( mp, node ) & 00229 * exp_igr % ez ( np, node ) ) ) 00230 fg_coulomb ( :, node ) = fg_coulomb ( :, node ) & 00231 + charge ( node ) * gauss * e_igdotr * pw_grid % g ( :, gpt ) 00232 END DO 00233 00234 ! compute the virial P*V 00235 denom = 1.0_dp / four_alpha_sq + 1.0_dp / pw_grid % gsq ( gpt ) 00236 IF (use_virial) THEN 00237 pv_g(1,1) = pv_g(1,1) + factor*(1.0_dp - 2.0_dp*pw_grid%g(1,gpt)*pw_grid%g(1,gpt)*denom) 00238 pv_g(1,2) = pv_g(1,2) - factor*( 2.0_dp*pw_grid%g(1,gpt)*pw_grid%g(2,gpt)*denom) 00239 pv_g(1,3) = pv_g(1,3) - factor*( 2.0_dp*pw_grid%g(1,gpt)*pw_grid%g(3,gpt)*denom) 00240 pv_g(2,1) = pv_g(2,1) - factor*( 2.0_dp*pw_grid%g(2,gpt)*pw_grid%g(1,gpt)*denom) 00241 pv_g(2,2) = pv_g(2,2) + factor*(1.0_dp - 2.0_dp*pw_grid%g(2,gpt)*pw_grid%g(2,gpt)*denom) 00242 pv_g(2,3) = pv_g(2,3) - factor*( 2.0_dp*pw_grid%g(2,gpt)*pw_grid%g(3,gpt)*denom) 00243 pv_g(3,1) = pv_g(3,1) - factor*( 2.0_dp*pw_grid%g(3,gpt)*pw_grid%g(1,gpt)*denom) 00244 pv_g(3,2) = pv_g(3,2) - factor*( 2.0_dp*pw_grid%g(3,gpt)*pw_grid%g(2,gpt)*denom) 00245 pv_g(3,3) = pv_g(3,3) + factor*(1.0_dp - 2.0_dp*pw_grid%g(3,gpt)*pw_grid%g(3,gpt)*denom) 00246 END IF 00247 END DO 00248 00249 vg_coulomb = vg_coulomb * pref 00250 IF (use_virial) pv_g = pv_g * pref 00251 IF (atenergy) e_coulomb = e_coulomb * pref 00252 00253 fg_coulomb = fg_coulomb * (2.0_dp*pref) 00254 00255 CALL structure_factor_deallocate ( exp_igr ) 00256 00257 DEALLOCATE ( charge, summe, STAT = isos ) 00258 CPPostcondition(isos==0,cp_failure_level,routineP,error,failure) 00259 NULLIFY ( charge ) 00260 00261 CALL timestop(handle) 00262 00263 END SUBROUTINE ewald_evaluate 00264 00265 ! ***************************************************************************** 00272 SUBROUTINE ewald_self ( ewald_env, cell, atomic_kind_set, local_particles, e_self,& 00273 e_neut, charges, error) 00274 00275 TYPE(ewald_environment_type), POINTER :: ewald_env 00276 TYPE(cell_type), POINTER :: cell 00277 TYPE(atomic_kind_type), POINTER :: atomic_kind_set( : ) 00278 TYPE(distribution_1d_type), POINTER :: local_particles 00279 REAL(KIND=dp), INTENT(OUT) :: e_self, e_neut 00280 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 00281 TYPE(cp_error_type), INTENT(inout) :: error 00282 00283 CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_self', 00284 routineP = moduleN//':'//routineN 00285 00286 INTEGER :: ewald_type, group, 00287 iparticle_kind, 00288 nparticle_kind, 00289 nparticle_local 00290 LOGICAL :: is_shell 00291 REAL(KIND=dp) :: alpha, mm_radius, q, q_neutg, 00292 q_self, q_sum, qcore, qshell 00293 TYPE(atomic_kind_type), POINTER :: atomic_kind 00294 TYPE(shell_kind_type), POINTER :: shell 00295 00296 CALL ewald_env_get ( ewald_env, ewald_type = ewald_type, & 00297 alpha = alpha, group = group ,error=error) 00298 q_neutg = 0.0_dp 00299 q_self = 0.0_dp 00300 q_sum = 0.0_dp 00301 nparticle_kind = SIZE ( atomic_kind_set ) 00302 IF (ASSOCIATED(charges)) THEN 00303 q_self = DOT_PRODUCT(charges,charges) 00304 q_sum = SUM(charges) 00305 ! check and abort.. 00306 DO iparticle_kind=1,nparticle_kind 00307 atomic_kind => atomic_kind_set(iparticle_kind) 00308 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius) 00309 IF(mm_radius>0.0_dp)THEN 00310 CALL cp_unimplemented_error(fromWhere=routineP, & 00311 message="Array of charges not implemented for mm_radius>0.0 !!",& 00312 error=error, error_level=cp_failure_level) 00313 END IF 00314 END DO 00315 ELSE 00316 DO iparticle_kind=1,nparticle_kind 00317 atomic_kind => atomic_kind_set(iparticle_kind) 00318 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius, & 00319 qeff=q, shell_active=is_shell, shell=shell) 00320 nparticle_local = local_particles%n_el(iparticle_kind) 00321 IF (is_shell) THEN 00322 CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell, error=error) 00323 ! MI: the core-shell ES interaction, when core and shell belong to the same ion, is excluded 00324 ! in the nonbond correction term. Therefore, here the self interaction is computed entirely 00325 q_self = q_self + qcore*qcore*nparticle_local + qshell*qshell*nparticle_local 00326 q_sum = q_sum + qcore*nparticle_local + qshell*nparticle_local 00327 IF (mm_radius > 0) THEN 00328 ! the core is always a point charge 00329 q_neutg = q_neutg + 2.0_dp*qshell*mm_radius**2 00330 END IF 00331 ELSE 00332 q_self = q_self + q*q*nparticle_local 00333 q_sum = q_sum + q*nparticle_local 00334 IF (mm_radius > 0) THEN 00335 q_neutg = q_neutg + 2.0_dp*q*mm_radius**2 00336 END IF 00337 END IF 00338 END DO 00339 00340 CALL mp_sum ( q_self, group ) 00341 CALL mp_sum ( q_sum, group ) 00342 END IF 00343 00344 e_neut = 0.0_dp 00345 e_self = 0.0_dp 00346 IF ( ewald_type /= do_ewald_none) THEN 00347 e_self = -q_self * alpha * oorootpi 00348 e_neut = -q_sum*pi/(2.0_dp*cell%deth)*(q_sum/alpha**2 - q_neutg) 00349 END IF 00350 00351 END SUBROUTINE ewald_self 00352 00353 ! ***************************************************************************** 00359 SUBROUTINE ewald_self_atom ( ewald_env, cell, atomic_kind_set, local_particles, e_self,& 00360 charges, error) 00361 00362 TYPE(ewald_environment_type), POINTER :: ewald_env 00363 TYPE(cell_type), POINTER :: cell 00364 TYPE(atomic_kind_type), DIMENSION(:), 00365 POINTER :: atomic_kind_set( : ) 00366 TYPE(distribution_1d_type), POINTER :: local_particles 00367 REAL(KIND=dp), DIMENSION(:), 00368 INTENT(INOUT) :: e_self(:) 00369 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 00370 TYPE(cp_error_type), INTENT(inout) :: error 00371 00372 CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_self_atom', 00373 routineP = moduleN//':'//routineN 00374 00375 INTEGER :: ewald_type, ii, 00376 iparticle_kind, 00377 iparticle_local, 00378 nparticle_kind, 00379 nparticle_local 00380 LOGICAL :: is_shell 00381 REAL(KIND=dp) :: alpha, fself, q, qcore, qshell 00382 TYPE(atomic_kind_type), POINTER :: atomic_kind 00383 TYPE(shell_kind_type), POINTER :: shell 00384 00385 CALL ewald_env_get(ewald_env,ewald_type=ewald_type,alpha=alpha,error=error) 00386 00387 fself = alpha*oorootpi 00388 00389 IF (ewald_type /= do_ewald_none) THEN 00390 nparticle_kind = SIZE(atomic_kind_set) 00391 IF (ASSOCIATED(charges)) THEN 00392 CALL cp_unimplemented_error(fromWhere=routineP, & 00393 message="Atomic energy not implemented for charges",& 00394 error=error, error_level=cp_failure_level) 00395 ELSE 00396 DO iparticle_kind=1,nparticle_kind 00397 atomic_kind => atomic_kind_set(iparticle_kind) 00398 nparticle_local = local_particles%n_el(iparticle_kind) 00399 CALL get_atomic_kind(atomic_kind=atomic_kind,qeff=q,& 00400 shell_active=is_shell,shell=shell) 00401 IF (is_shell) THEN 00402 CALL get_shell(shell=shell,charge_core=qcore,charge_shell=qshell,error=error) 00403 DO iparticle_local=1,nparticle_local 00404 ii = local_particles%list(iparticle_kind)%array(iparticle_local) 00405 e_self(ii) = e_self(ii) - (qcore*qcore + qshell*qshell)*fself 00406 END DO 00407 ELSE 00408 DO iparticle_local=1,nparticle_local 00409 ii = local_particles%list(iparticle_kind)%array(iparticle_local) 00410 e_self(ii) = e_self(ii) - q*q*fself 00411 END DO 00412 END IF 00413 END DO 00414 END IF 00415 END IF 00416 00417 END SUBROUTINE ewald_self_atom 00418 00419 ! ***************************************************************************** 00424 SUBROUTINE ewald_print ( iw, pot_nonbond, e_gspace, e_self, e_neut, e_bonded ) 00425 00426 INTEGER, INTENT(IN) :: iw 00427 REAL(KIND=dp), INTENT(IN) :: pot_nonbond, e_gspace, 00428 e_self, e_neut, e_bonded 00429 00430 CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_print', 00431 routineP = moduleN//':'//routineN 00432 00433 IF (iw>0) THEN 00434 WRITE ( iw, '( A, A )' ) ' *********************************', & 00435 '**********************************************' 00436 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' INITIAL GSPACE ENERGY', & 00437 '[hartree]', '= ', e_gspace 00438 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' INITIAL NONBONDED ENERGY', & 00439 '[hartree]', '= ', pot_nonbond 00440 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' SELF ENERGY CORRECTION', & 00441 '[hartree]', '= ', e_self 00442 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' NEUT. BACKGROUND',& 00443 '[hartree]', '= ', e_neut 00444 WRITE ( iw, '( A, A, T35, A, T56, E25.15 )' ) ' BONDED CORRECTION',& 00445 '[hartree]','= ', e_bonded 00446 WRITE ( iw, '( A, A )' ) ' *********************************', & 00447 '**********************************************' 00448 END IF 00449 END SUBROUTINE ewald_print 00450 00451 END MODULE ewalds
1.7.3