|
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 SUBROUTINE debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, & 00012 particle_set, local_particles, iw, debug_r_space, error ) 00013 USE cell_types, ONLY: cell_type 00014 USE distribution_1d_types, ONLY: distribution_1d_type 00015 USE ewald_environment_types, ONLY: ewald_environment_type 00016 USE ewald_pw_types, ONLY: ewald_pw_type 00017 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 00018 neighbor_kind_pairs_type 00019 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 00020 fist_nonbond_env_type,& 00021 pos_type 00022 USE kinds, ONLY: dp 00023 USE mathlib, ONLY: matvec_3x3 00024 USE parallel_rng_types, ONLY: UNIFORM,& 00025 create_rng_stream,& 00026 delete_rng_stream,& 00027 next_random_number,& 00028 rng_stream_type,& 00029 random_numbers 00030 USE particle_types, ONLY: particle_type 00031 USE ewalds_multipole, ONLY: ewald_multipole_evaluate 00032 #include "cp_common_uses.h" 00033 IMPLICIT NONE 00034 00035 TYPE charge_mono_type 00036 REAL(KIND=dp), DIMENSION(:), 00037 POINTER :: charge 00038 REAL(KIND=dp), DIMENSION(:,:), 00039 POINTER :: pos 00040 END TYPE charge_mono_type 00041 TYPE multi_charge_type 00042 TYPE(charge_mono_type), DIMENSION(:), 00043 POINTER :: charge_typ 00044 END TYPE multi_charge_type 00045 TYPE(ewald_environment_type), POINTER :: ewald_env 00046 TYPE(ewald_pw_type), POINTER :: ewald_pw 00047 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00048 TYPE(cell_type), POINTER :: cell 00049 TYPE(particle_type), DIMENSION(:), 00050 POINTER :: particle_set 00051 TYPE(distribution_1d_type), POINTER :: local_particles 00052 INTEGER, INTENT(IN) :: iw 00053 LOGICAL, INTENT(IN) :: debug_r_space 00054 TYPE(cp_error_type), INTENT(inout) :: error 00055 00056 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipoles', 00057 routineP = "ewalds_multipole_debug"//':'//routineN 00058 00059 INTEGER :: i, ind, nparticles, stat 00060 LOGICAL :: failure 00061 LOGICAL, DIMENSION(3) :: task 00062 REAL(KIND=dp) :: e_neut, e_self, g_energy, 00063 r_energy, debug_energy 00064 REAL(KIND=dp), POINTER, DIMENSION(:) :: charges 00065 REAL(KIND=dp), POINTER, 00066 DIMENSION(:, :) :: dipoles, g_forces, g_pv, 00067 r_forces, r_pv, e_field1, 00068 e_field2 00069 REAL(KIND=dp), POINTER, 00070 DIMENSION(:, :, :) :: quadrupoles 00071 TYPE(rng_stream_type), POINTER :: random_stream 00072 TYPE(multi_charge_type), DIMENSION(:), 00073 POINTER :: multipoles 00074 00075 failure = .FALSE. 00076 NULLIFY(random_stream, multipoles, charges, dipoles, g_forces, g_pv,& 00077 r_forces, r_pv, e_field1, e_field2) 00078 CALL create_rng_stream(random_stream,name="DEBUG_EWALD_MULTIPOLE",& 00079 distribution_type=UNIFORM,error=error) 00080 ! check: charge - charge 00081 task = .FALSE. 00082 nparticles = SIZE(particle_set) 00083 00084 ! Allocate charges, dipoles, quadrupoles 00085 ALLOCATE(charges(nparticles),stat=stat) 00086 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00087 ALLOCATE(dipoles(3,nparticles),stat=stat) 00088 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00089 ALLOCATE(quadrupoles(3,3,nparticles),stat=stat) 00090 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00091 00092 ! Allocate arrays for forces 00093 ALLOCATE(r_forces(3,nparticles),stat=stat) 00094 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00095 ALLOCATE(g_forces(3,nparticles),stat=stat) 00096 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00097 ALLOCATE(e_field1(3,nparticles),stat=stat) 00098 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00099 ALLOCATE(e_field2(3,nparticles),stat=stat) 00100 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00101 ALLOCATE(g_pv(3,3),stat=stat) 00102 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00103 ALLOCATE(r_pv(3,3),stat=stat) 00104 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00105 00106 ! Debug CHARGES-CHARGES 00107 task(1) = .TRUE. 00108 charges = 0.0_dp 00109 dipoles = 0.0_dp 00110 quadrupoles = 0.0_dp 00111 r_forces = 0.0_dp 00112 g_forces = 0.0_dp 00113 e_field1 = 0.0_dp 00114 e_field2 = 0.0_dp 00115 g_pv = 0.0_dp 00116 r_pv = 0.0_dp 00117 g_energy = 0.0_dp 00118 r_energy = 0.0_dp 00119 e_neut = 0.0_dp 00120 e_self = 0.0_dp 00121 00122 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, & 00123 random_stream=random_stream, particle_set=particle_set, charges=charges,& 00124 error=error) 00125 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "CHARGE", echarge=1.0_dp, & 00126 random_stream=random_stream, particle_set=particle_set, charges=charges,& 00127 error=error) 00128 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy,& 00129 debug_r_space, error) 00130 00131 WRITE(*,*)"DEBUG ENERGY (CHARGE-CHARGE): ", debug_energy 00132 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 00133 particle_set, local_particles, g_energy, r_energy, e_neut, e_self,& 00134 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE.,& 00135 charges=charges,dipoles=dipoles,quadrupoles=quadrupoles,forces_local=g_forces,& 00136 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.,& 00137 error=error) 00138 CALL release_multi_type(multipoles, error=error) 00139 00140 00141 ! Debug CHARGES-DIPOLES 00142 task(1) = .TRUE. 00143 task(2) = .TRUE. 00144 charges = 0.0_dp 00145 dipoles = 0.0_dp 00146 quadrupoles = 0.0_dp 00147 r_forces = 0.0_dp 00148 g_forces = 0.0_dp 00149 e_field1 = 0.0_dp 00150 e_field2 = 0.0_dp 00151 g_pv = 0.0_dp 00152 r_pv = 0.0_dp 00153 g_energy = 0.0_dp 00154 r_energy = 0.0_dp 00155 e_neut = 0.0_dp 00156 e_self = 0.0_dp 00157 00158 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, & 00159 random_stream=random_stream, particle_set=particle_set, charges=charges,& 00160 error=error) 00161 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "DIPOLE", echarge=0.5_dp, & 00162 random_stream=random_stream, particle_set=particle_set, dipoles=dipoles,& 00163 error=error) 00164 WRITE(*,'("CHARGES",F15.9)')charges 00165 WRITE(*,'("DIPOLES",3F15.9)')dipoles 00166 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 00167 debug_r_space, error) 00168 00169 WRITE(*,*)"DEBUG ENERGY (CHARGE-DIPOLE): ", debug_energy 00170 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 00171 particle_set, local_particles, g_energy, r_energy, e_neut, e_self,& 00172 task,do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE.,& 00173 charges=charges,dipoles=dipoles,quadrupoles=quadrupoles,forces_local=g_forces,& 00174 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.,& 00175 error=error) 00176 CALL release_multi_type(multipoles, error=error) 00177 00178 ! Debug DIPOLES-DIPOLES 00179 task(2) = .TRUE. 00180 charges = 0.0_dp 00181 dipoles = 0.0_dp 00182 quadrupoles = 0.0_dp 00183 r_forces = 0.0_dp 00184 g_forces = 0.0_dp 00185 e_field1 = 0.0_dp 00186 e_field2 = 0.0_dp 00187 g_pv = 0.0_dp 00188 r_pv = 0.0_dp 00189 g_energy = 0.0_dp 00190 r_energy = 0.0_dp 00191 e_neut = 0.0_dp 00192 e_self = 0.0_dp 00193 00194 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "DIPOLE", echarge=10000.0_dp, & 00195 random_stream=random_stream, particle_set=particle_set, dipoles=dipoles,& 00196 error=error) 00197 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "DIPOLE", echarge=20000._dp, & 00198 random_stream=random_stream, particle_set=particle_set, dipoles=dipoles,& 00199 error=error) 00200 WRITE(*,'("DIPOLES",3F15.9)')dipoles 00201 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 00202 debug_r_space, error) 00203 00204 WRITE(*,*)"DEBUG ENERGY (DIPOLE-DIPOLE): ", debug_energy 00205 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 00206 particle_set, local_particles, g_energy, r_energy, e_neut, e_self,& 00207 task,do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE.,& 00208 charges=charges,dipoles=dipoles,quadrupoles=quadrupoles,forces_local=g_forces,& 00209 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.,& 00210 error=error) 00211 CALL release_multi_type(multipoles, error=error) 00212 00213 ! Debug CHARGES-QUADRUPOLES 00214 task(1) = .TRUE. 00215 task(3) = .TRUE. 00216 charges = 0.0_dp 00217 dipoles = 0.0_dp 00218 quadrupoles = 0.0_dp 00219 r_forces = 0.0_dp 00220 g_forces = 0.0_dp 00221 e_field1 = 0.0_dp 00222 e_field2 = 0.0_dp 00223 g_pv = 0.0_dp 00224 r_pv = 0.0_dp 00225 g_energy = 0.0_dp 00226 r_energy = 0.0_dp 00227 e_neut = 0.0_dp 00228 e_self = 0.0_dp 00229 00230 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, & 00231 random_stream=random_stream, particle_set=particle_set, charges=charges,& 00232 error=error) 00233 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "QUADRUPOLE", echarge=10.0_dp, & 00234 random_stream=random_stream, particle_set=particle_set, quadrupoles=quadrupoles,& 00235 error=error) 00236 WRITE(*,'("CHARGES",F15.9)')charges 00237 WRITE(*,'("QUADRUPOLES",9F15.9)')quadrupoles 00238 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 00239 debug_r_space, error) 00240 00241 WRITE(*,*)"DEBUG ENERGY (CHARGE-QUADRUPOLE): ", debug_energy 00242 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 00243 particle_set, local_particles, g_energy, r_energy, e_neut, e_self,& 00244 task,do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE.,& 00245 charges=charges,dipoles=dipoles,quadrupoles=quadrupoles,forces_local=g_forces,& 00246 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.,& 00247 error=error) 00248 CALL release_multi_type(multipoles, error=error) 00249 00250 ! Debug DIPOLES-QUADRUPOLES 00251 task(2) = .TRUE. 00252 task(3) = .TRUE. 00253 charges = 0.0_dp 00254 dipoles = 0.0_dp 00255 quadrupoles = 0.0_dp 00256 r_forces = 0.0_dp 00257 g_forces = 0.0_dp 00258 e_field1 = 0.0_dp 00259 e_field2 = 0.0_dp 00260 g_pv = 0.0_dp 00261 r_pv = 0.0_dp 00262 g_energy = 0.0_dp 00263 r_energy = 0.0_dp 00264 e_neut = 0.0_dp 00265 e_self = 0.0_dp 00266 00267 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "DIPOLE", echarge=10000.0_dp, & 00268 random_stream=random_stream, particle_set=particle_set, dipoles=dipoles,& 00269 error=error) 00270 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "QUADRUPOLE", echarge=10000.0_dp, & 00271 random_stream=random_stream, particle_set=particle_set, quadrupoles=quadrupoles,& 00272 error=error) 00273 WRITE(*,'("DIPOLES",3F15.9)')dipoles 00274 WRITE(*,'("QUADRUPOLES",9F15.9)')quadrupoles 00275 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 00276 debug_r_space, error) 00277 00278 WRITE(*,*)"DEBUG ENERGY (DIPOLE-QUADRUPOLE): ", debug_energy 00279 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 00280 particle_set, local_particles, g_energy, r_energy, e_neut, e_self,& 00281 task,do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE.,& 00282 charges=charges,dipoles=dipoles,quadrupoles=quadrupoles,forces_local=g_forces,& 00283 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.,& 00284 error=error) 00285 CALL release_multi_type(multipoles, error=error) 00286 00287 ! Debug QUADRUPOLES-QUADRUPOLES 00288 task(3) = .TRUE. 00289 charges = 0.0_dp 00290 dipoles = 0.0_dp 00291 quadrupoles = 0.0_dp 00292 r_forces = 0.0_dp 00293 g_forces = 0.0_dp 00294 e_field1 = 0.0_dp 00295 e_field2 = 0.0_dp 00296 g_pv = 0.0_dp 00297 r_pv = 0.0_dp 00298 g_energy = 0.0_dp 00299 r_energy = 0.0_dp 00300 e_neut = 0.0_dp 00301 e_self = 0.0_dp 00302 00303 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "QUADRUPOLE", echarge=-20000.0_dp, & 00304 random_stream=random_stream, particle_set=particle_set, quadrupoles=quadrupoles,& 00305 error=error) 00306 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "QUADRUPOLE", echarge=10000.0_dp, & 00307 random_stream=random_stream, particle_set=particle_set, quadrupoles=quadrupoles,& 00308 error=error) 00309 WRITE(*,'("QUADRUPOLES",9F15.9)')quadrupoles 00310 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 00311 debug_r_space, error) 00312 00313 WRITE(*,*)"DEBUG ENERGY (QUADRUPOLE-QUADRUPOLE): ", debug_energy 00314 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 00315 particle_set, local_particles, g_energy, r_energy, e_neut, e_self,& 00316 task,do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE.,& 00317 charges=charges,dipoles=dipoles,quadrupoles=quadrupoles,forces_local=g_forces,& 00318 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.,& 00319 error=error) 00320 CALL release_multi_type(multipoles, error=error) 00321 00322 00323 DEALLOCATE(charges,stat=stat) 00324 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00325 DEALLOCATE(dipoles,stat=stat) 00326 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00327 DEALLOCATE(quadrupoles,stat=stat) 00328 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00329 DEALLOCATE(r_forces,stat=stat) 00330 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00331 DEALLOCATE(g_forces,stat=stat) 00332 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00333 DEALLOCATE(e_field1,stat=stat) 00334 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00335 DEALLOCATE(e_field2,stat=stat) 00336 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00337 DEALLOCATE(g_pv,stat=stat) 00338 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00339 DEALLOCATE(r_pv,stat=stat) 00340 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00341 CALL delete_rng_stream(random_stream,error=error) 00342 00343 CONTAINS 00344 ! ***************************************************************************** 00349 SUBROUTINE debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles,& 00350 energy, debug_r_space, error) 00351 IMPLICIT NONE 00352 TYPE(particle_type), DIMENSION(:), 00353 POINTER :: particle_set 00354 TYPE(cell_type), POINTER :: cell 00355 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00356 TYPE(multi_charge_type), DIMENSION(:), 00357 POINTER :: multipoles 00358 REAL(KIND=dp), INTENT(OUT) :: energy 00359 LOGICAL, INTENT(IN) :: debug_r_space 00360 TYPE(cp_error_type), INTENT(inout) :: error 00361 00362 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipole_low', 00363 routineP = "ewalds_multipole_debug"//':'//routineN 00364 00365 LOGICAL :: failure 00366 INTEGER :: atom_a, atom_b, iend, igrp, 00367 ikind, ilist, ipair, istart, 00368 jkind, nkinds, npairs,k,k1,l,l1, 00369 ncells, icell, jcell, kcell 00370 INTEGER, DIMENSION(:, :), POINTER :: list 00371 REAL(KIND=dp) :: r, rab2, rab2_max, fac_ij, q 00372 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi, rab, rm, rab0 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 energy = 0.0_dp 00380 CALL fist_nonbond_env_get (nonbond_env, nonbonded=nonbonded, natom_types = nkinds,& 00381 r_last_update=r_last_update,r_last_update_pbc=r_last_update_pbc, error=error) 00382 rab2_max = HUGE(0.0_dp) 00383 IF (debug_r_space) THEN 00384 ! This debugs the real space part of the multipole Ewald summation scheme 00385 ! Starting the force loop 00386 Lists: DO ilist=1,nonbonded%nlists 00387 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 00388 npairs=neighbor_kind_pair%npairs 00389 IF (npairs ==0) CYCLE 00390 list => neighbor_kind_pair%list 00391 cvi = neighbor_kind_pair%cell_vector 00392 CALL matvec_3x3(cell_v, cell%hmat, cvi) 00393 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 00394 istart = neighbor_kind_pair%grp_kind_start(igrp) 00395 iend = neighbor_kind_pair%grp_kind_end(igrp) 00396 ikind = neighbor_kind_pair%ij_kind(1,igrp) 00397 jkind = neighbor_kind_pair%ij_kind(2,igrp) 00398 Pairs: DO ipair = istart, iend 00399 fac_ij = 1.0_dp 00400 atom_a = list(1,ipair) 00401 atom_b = list(2,ipair) 00402 IF (atom_a==atom_b) fac_ij = 0.5_dp 00403 rab = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 00404 rab = rab + cell_v 00405 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2 00406 IF (rab2 <= rab2_max) THEN 00407 00408 DO k = 1, SIZE(multipoles(atom_a)%charge_typ) 00409 DO k1 = 1, SIZE(multipoles(atom_a)%charge_typ(k)%charge) 00410 00411 DO l = 1, SIZE(multipoles(atom_b)%charge_typ) 00412 DO l1 = 1, SIZE(multipoles(atom_b)%charge_typ(l)%charge) 00413 00414 rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:,l1) - multipoles(atom_a)%charge_typ(k)%pos(:,k1) 00415 r = SQRT(DOT_PRODUCT(rm,rm)) 00416 q = multipoles(atom_b)%charge_typ(l)%charge(l1) * multipoles(atom_a)%charge_typ(k)%charge(k1) 00417 energy = energy + q / r * fac_ij 00418 END DO 00419 END DO 00420 00421 END DO 00422 END DO 00423 00424 END IF 00425 END DO Pairs 00426 END DO Kind_Group_Loop 00427 END DO Lists 00428 ELSE 00429 ncells = 6 00430 !Debugs the sum of real + space terms.. (Charge-Charge and Charge-Dipole should be anyway wrong but 00431 !all the other terms should be correct) 00432 DO atom_a = 1, SIZE(particle_set) 00433 DO atom_b = atom_a, SIZE(particle_set) 00434 fac_ij = 1.0_dp 00435 IF (atom_a==atom_b) fac_ij = 0.5_dp 00436 rab0 = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 00437 ! Loop over cells 00438 DO icell = -ncells,ncells 00439 DO jcell = -ncells,ncells 00440 DO kcell = -ncells,ncells 00441 cell_v= MATMUL(cell%hmat,REAL((/icell,jcell,kcell/),KIND=dp)) 00442 IF (ALL(cell_v==0.0_dp).AND.(atom_a==atom_b)) CYCLE 00443 rab = rab0 + cell_v 00444 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2 00445 IF (rab2 <= rab2_max) THEN 00446 00447 DO k = 1, SIZE(multipoles(atom_a)%charge_typ) 00448 DO k1 = 1, SIZE(multipoles(atom_a)%charge_typ(k)%charge) 00449 00450 DO l = 1, SIZE(multipoles(atom_b)%charge_typ) 00451 DO l1 = 1, SIZE(multipoles(atom_b)%charge_typ(l)%charge) 00452 00453 rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:,l1) - multipoles(atom_a)%charge_typ(k)%pos(:,k1) 00454 r = SQRT(DOT_PRODUCT(rm,rm)) 00455 q = multipoles(atom_b)%charge_typ(l)%charge(l1) * multipoles(atom_a)%charge_typ(k)%charge(k1) 00456 energy = energy + q / r * fac_ij 00457 END DO 00458 END DO 00459 00460 END DO 00461 END DO 00462 00463 END IF 00464 END DO 00465 END DO 00466 END DO 00467 END DO 00468 END DO 00469 END IF 00470 END SUBROUTINE debug_ewald_multipole_low 00471 00472 ! ***************************************************************************** 00477 SUBROUTINE create_multi_type(multipoles, idim, istart, iend, label, echarge,& 00478 random_stream, particle_set, charges, dipoles, quadrupoles, error) 00479 IMPLICIT NONE 00480 TYPE(multi_charge_type), DIMENSION(:), 00481 POINTER :: multipoles 00482 INTEGER, INTENT(IN) :: idim, istart, iend 00483 CHARACTER(LEN=*), INTENT(IN) :: label 00484 REAL(KIND=dp), INTENT(IN) :: echarge 00485 TYPE(rng_stream_type), POINTER :: random_stream 00486 TYPE(particle_type), DIMENSION(:), 00487 POINTER :: particle_set 00488 REAL(KIND=dp), DIMENSION(:), POINTER, 00489 OPTIONAL :: charges 00490 REAL(KIND=dp), DIMENSION(:,:), POINTER, 00491 OPTIONAL :: dipoles 00492 REAL(KIND=dp), DIMENSION(:,:,:), POINTER, 00493 OPTIONAL :: quadrupoles 00494 TYPE(cp_error_type), INTENT(inout) :: error 00495 00496 CHARACTER(len=*), PARAMETER :: routineN = 'create_multi_type', 00497 routineP = "ewalds_multipole_debug"//':'//routineN 00498 00499 LOGICAL :: failure 00500 INTEGER :: i, stat, isize, k, l, m 00501 REAL(KIND=dp) :: dx, rvec(3), rvec1(3), rvec2(3), r2 00502 00503 failure = .FALSE. 00504 IF (ASSOCIATED(multipoles)) THEN 00505 CPPostcondition(SIZE(multipoles)==idim,cp_failure_level,routineP,error,failure) 00506 ELSE 00507 ALLOCATE(multipoles(idim),stat=stat) 00508 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00509 DO i = 1, idim 00510 NULLIFY(multipoles(i)%charge_typ) 00511 END DO 00512 END IF 00513 DO i = istart, iend 00514 IF (ASSOCIATED(multipoles(i)%charge_typ)) THEN 00515 ! make a copy of the array and enlarge the array type by 1 00516 isize = SIZE(multipoles(i)%charge_typ)+1 00517 ELSE 00518 isize = 1 00519 END IF 00520 CALL reallocate_charge_type(multipoles(i)%charge_typ,1,isize,error) 00521 SELECT CASE(label) 00522 CASE("CHARGE") 00523 CPPostcondition(PRESENT(charges),cp_failure_level,routineP,error,failure) 00524 CPPostcondition(ASSOCIATED(charges),cp_failure_level,routineP,error,failure) 00525 ALLOCATE(multipoles(i)%charge_typ(isize)%charge(1),stat=stat) 00526 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00527 ALLOCATE(multipoles(i)%charge_typ(isize)%pos(3,1),stat=stat) 00528 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00529 00530 multipoles(i)%charge_typ(isize)%charge(1) = echarge 00531 multipoles(i)%charge_typ(isize)%pos(1:3,1) = 0.0_dp 00532 charges(i) = charges(i) + echarge 00533 CASE("DIPOLE") 00534 dx = 1.0E-4_dp 00535 CPPostcondition(PRESENT(dipoles),cp_failure_level,routineP,error,failure) 00536 CPPostcondition(ASSOCIATED(dipoles),cp_failure_level,routineP,error,failure) 00537 ALLOCATE(multipoles(i)%charge_typ(isize)%charge(2),stat=stat) 00538 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00539 ALLOCATE(multipoles(i)%charge_typ(isize)%pos(3,2),stat=stat) 00540 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00541 CALL random_numbers(rvec, random_stream, error) 00542 rvec = rvec/(2.0_dp*SQRT(DOT_PRODUCT(rvec,rvec)))*dx 00543 multipoles(i)%charge_typ(isize)%charge(1) = echarge 00544 multipoles(i)%charge_typ(isize)%pos(1:3,1) = rvec 00545 multipoles(i)%charge_typ(isize)%charge(2) =-echarge 00546 multipoles(i)%charge_typ(isize)%pos(1:3,2) =-rvec 00547 00548 dipoles(:,i) = dipoles(:,i) + 2.0_dp*echarge*rvec 00549 CASE("QUADRUPOLE") 00550 dx = 1.0E-2_dp 00551 CPPostcondition(PRESENT(quadrupoles),cp_failure_level,routineP,error,failure) 00552 CPPostcondition(ASSOCIATED(quadrupoles),cp_failure_level,routineP,error,failure) 00553 ALLOCATE(multipoles(i)%charge_typ(isize)%charge(4),stat=stat) 00554 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00555 ALLOCATE(multipoles(i)%charge_typ(isize)%pos(3,4),stat=stat) 00556 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00557 CALL random_numbers(rvec1, random_stream, error) 00558 CALL random_numbers(rvec2, random_stream, error) 00559 rvec1 = rvec1/SQRT(DOT_PRODUCT(rvec1,rvec1)) 00560 rvec2 = rvec2 - DOT_PRODUCT(rvec2,rvec1)*rvec1 00561 rvec2 = rvec2/SQRT(DOT_PRODUCT(rvec2,rvec2)) 00562 ! 00563 rvec1 = rvec1/2.0_dp*dx 00564 rvec2 = rvec2/2.0_dp*dx 00565 ! + (4) ^ - (1) 00566 ! |rvec2 00567 ! | 00568 ! 0------> rvec1 00569 ! 00570 ! 00571 ! - (3) + (2) 00572 multipoles(i)%charge_typ(isize)%charge(1) = -echarge 00573 multipoles(i)%charge_typ(isize)%pos(1:3,1) = rvec1+rvec2 00574 multipoles(i)%charge_typ(isize)%charge(2) = echarge 00575 multipoles(i)%charge_typ(isize)%pos(1:3,2) = rvec1-rvec2 00576 multipoles(i)%charge_typ(isize)%charge(3) = -echarge 00577 multipoles(i)%charge_typ(isize)%pos(1:3,3) = -rvec1-rvec2 00578 multipoles(i)%charge_typ(isize)%charge(4) = echarge 00579 multipoles(i)%charge_typ(isize)%pos(1:3,4) = -rvec1+rvec2 00580 00581 DO k = 1, 4 00582 r2 = DOT_PRODUCT(multipoles(i)%charge_typ(isize)%pos(:,k),multipoles(i)%charge_typ(isize)%pos(:,k)) 00583 DO l = 1, 3 00584 DO m = 1, 3 00585 quadrupoles(m,l,i) = quadrupoles(m,l,i) + 3.0_dp* 0.5_dp *multipoles(i)%charge_typ(isize)%charge(k) * & 00586 multipoles(i)%charge_typ(isize)%pos(l,k) * & 00587 multipoles(i)%charge_typ(isize)%pos(m,k) 00588 IF (m==l) quadrupoles(m,l,i) = quadrupoles(m,l,i) - 0.5_dp * multipoles(i)%charge_typ(isize)%charge(k) * r2 00589 END DO 00590 END DO 00591 END DO 00592 00593 END SELECT 00594 END DO 00595 END SUBROUTINE create_multi_type 00596 00597 ! ***************************************************************************** 00602 SUBROUTINE release_multi_type(multipoles, error) 00603 IMPLICIT NONE 00604 TYPE(multi_charge_type), DIMENSION(:), 00605 POINTER :: multipoles 00606 TYPE(cp_error_type), INTENT(inout) :: error 00607 00608 CHARACTER(len=*), PARAMETER :: routineN = 'release_multi_type', 00609 routineP = "ewalds_multipole_debug"//':'//routineN 00610 00611 LOGICAL :: failure 00612 INTEGER :: i, stat, j 00613 00614 failure = .FALSE. 00615 IF (ASSOCIATED(multipoles)) THEN 00616 DO i = 1, SIZE(multipoles) 00617 DO j = 1, SIZE(multipoles(i)%charge_typ) 00618 DEALLOCATE(multipoles(i)%charge_typ(j)%charge,stat=stat) 00619 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00620 DEALLOCATE(multipoles(i)%charge_typ(j)%pos,stat=stat) 00621 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00622 END DO 00623 DEALLOCATE(multipoles(i)%charge_typ,stat=stat) 00624 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00625 END DO 00626 END IF 00627 END SUBROUTINE release_multi_type 00628 00629 ! ***************************************************************************** 00634 SUBROUTINE reallocate_charge_type(charge_typ, istart, iend, error) 00635 IMPLICIT NONE 00636 TYPE(charge_mono_type), DIMENSION(:), 00637 POINTER :: charge_typ 00638 INTEGER, INTENT(IN) :: istart, iend 00639 TYPE(cp_error_type), INTENT(inout) :: error 00640 00641 CHARACTER(len=*), PARAMETER :: routineN = 'reallocate_charge_type', 00642 routineP = "ewalds_multipole_debug"//':'//routineN 00643 00644 LOGICAL :: failure 00645 INTEGER :: i, stat, isize, jsize, jsize1, jsize2, j 00646 TYPE(charge_mono_type), DIMENSION(:), 00647 POINTER :: charge_typ_bk 00648 00649 failure = .FALSE. 00650 IF (ASSOCIATED(charge_typ)) THEN 00651 isize = SIZE(charge_typ) 00652 ALLOCATE(charge_typ_bk(1:isize),stat=stat) 00653 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00654 DO j = 1, isize 00655 jsize = SIZE(charge_typ(j)%charge) 00656 ALLOCATE(charge_typ_bk(j)%charge(jsize),stat=stat) 00657 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00658 jsize1 = SIZE(charge_typ(j)%pos,1) 00659 jsize2 = SIZE(charge_typ(j)%pos,2) 00660 ALLOCATE(charge_typ_bk(j)%pos(jsize1,jsize2),stat=stat) 00661 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00662 charge_typ_bk(j)%pos = charge_typ(j)%pos 00663 charge_typ_bk(j)%charge = charge_typ(j)%charge 00664 END DO 00665 DO j = 1, SIZE(charge_typ) 00666 DEALLOCATE(charge_typ(j)%charge,stat=stat) 00667 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00668 DEALLOCATE(charge_typ(j)%pos,stat=stat) 00669 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00670 END DO 00671 DEALLOCATE(charge_typ,stat=stat) 00672 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00673 ! Reallocate 00674 ALLOCATE(charge_typ_bk(istart:iend),stat=stat) 00675 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00676 DO i = istart, isize 00677 jsize = SIZE(charge_typ_bk(j)%charge) 00678 ALLOCATE(charge_typ(j)%charge(jsize),stat=stat) 00679 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00680 jsize1 = SIZE(charge_typ_bk(j)%pos,1) 00681 jsize2 = SIZE(charge_typ_bk(j)%pos,2) 00682 ALLOCATE(charge_typ(j)%pos(jsize1,jsize2),stat=stat) 00683 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00684 charge_typ(j)%pos = charge_typ_bk(j)%pos 00685 charge_typ(j)%charge = charge_typ_bk(j)%charge 00686 END DO 00687 DO j = 1, SIZE(charge_typ_bk) 00688 DEALLOCATE(charge_typ_bk(j)%charge,stat=stat) 00689 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00690 DEALLOCATE(charge_typ_bk(j)%pos,stat=stat) 00691 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00692 END DO 00693 DEALLOCATE(charge_typ_bk,stat=stat) 00694 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00695 ELSE 00696 ALLOCATE(charge_typ(istart:iend), stat=stat) 00697 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00698 END IF 00699 00700 END SUBROUTINE reallocate_charge_type 00701 00702 END SUBROUTINE debug_ewald_multipoles 00703 00704 ! ***************************************************************************** 00709 SUBROUTINE debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, cell,& 00710 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw, & 00711 atomic_kind_set, mm_section, error ) 00712 USE atomic_kind_types, ONLY: atomic_kind_type 00713 USE cell_types, ONLY: cell_type 00714 USE distribution_1d_types, ONLY: distribution_1d_type 00715 USE ewald_environment_types, ONLY: ewald_environment_type 00716 USE ewald_pw_types, ONLY: ewald_pw_type 00717 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 00718 neighbor_kind_pairs_type 00719 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 00720 fist_nonbond_env_type,& 00721 pos_type 00722 USE input_section_types, ONLY: section_vals_type 00723 USE kinds, ONLY: dp 00724 USE particle_types, ONLY: particle_type 00725 USE ewalds_multipole, ONLY: ewald_multipole_evaluate 00726 USE fist_neighbor_list_control, ONLY: list_control 00727 #include "cp_common_uses.h" 00728 IMPLICIT NONE 00729 TYPE(ewald_environment_type), POINTER :: ewald_env 00730 TYPE(ewald_pw_type), POINTER :: ewald_pw 00731 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00732 TYPE(cell_type), POINTER :: cell 00733 TYPE(particle_type), POINTER :: particle_set(:) 00734 TYPE(distribution_1d_type), POINTER :: local_particles 00735 REAL(KIND=dp), DIMENSION(:), 00736 POINTER, OPTIONAL :: radii, charges 00737 REAL(KIND=dp), DIMENSION(:, :), 00738 POINTER, OPTIONAL :: dipoles 00739 REAL(KIND=dp), DIMENSION(:, :, :), 00740 POINTER, OPTIONAL :: quadrupoles 00741 LOGICAL, DIMENSION(3), INTENT(IN) :: task 00742 INTEGER, INTENT(IN) :: iw 00743 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:) 00744 TYPE(section_vals_type), POINTER :: mm_section 00745 TYPE(cp_error_type), INTENT(inout) :: error 00746 00747 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipoles_fields', 00748 routineP = "ewalds_multipole_debug"//':'//routineN 00749 00750 REAL(KIND=dp), POINTER, DIMENSION(:) :: lcharges, efield0 00751 REAL(KIND=dp), POINTER, DIMENSION(:,:) :: ldipoles 00752 REAL(KIND=dp), POINTER, DIMENSION(:,:) :: lquadrupoles 00753 REAL(KIND=dp) :: dq, dr, energy_local, energy_glob, e_neut, 00754 e_self, pv_local(3,3), pv_glob(3,3), ene(2), pot, coord(3), efield1n(3), tot_ene, 00755 o_tot_ene, enev(3,2), efield2n(3,3) 00756 REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE 00757 :: forces_local, forces_glob, efield1, efield2 00758 INTEGER :: nparticles,i,k,j,nparticle_local, ind, iparticle_kind 00759 TYPE(particle_type), POINTER, DIMENSION(:) :: shell_particle_set, core_particle_set 00760 TYPE(cp_logger_type), POINTER :: logger 00761 00762 NULLIFY(lcharges, ldipoles, lquadrupoles, shell_particle_set, core_particle_set) 00763 NULLIFY(logger) 00764 logger => cp_error_get_logger(error) 00765 00766 nparticles = SIZE(particle_set) 00767 nparticle_local = 0 00768 DO iparticle_kind=1,SIZE(local_particles%n_el) 00769 nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind) 00770 END DO 00771 ALLOCATE(lcharges(nparticles)) 00772 ALLOCATE(forces_glob(3,nparticles)) 00773 ALLOCATE(forces_local(3,nparticle_local)) 00774 ALLOCATE(efield0(nparticles)) 00775 ALLOCATE(efield1(3,nparticles)) 00776 ALLOCATE(efield2(9,nparticles)) 00777 forces_glob = 0.0_dp 00778 forces_local= 0.0_dp 00779 efield0 = 0.0_dp 00780 efield1 = 0.0_dp 00781 efield2 = 0.0_dp 00782 pv_local = 0.0_dp 00783 pv_glob = 0.0_dp 00784 energy_glob = 0.0_dp 00785 energy_local= 0.0_dp 00786 e_neut = 0.0_dp 00787 e_self = 0.0_dp 00788 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set,& 00789 local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .TRUE., .TRUE.,& 00790 .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob,& 00791 efield0, efield1, efield2, iw, do_debug=.FALSE., error=error ) 00792 o_tot_ene = energy_local + energy_glob + e_neut + e_self 00793 WRITE(*,*)"TOTAL ENERGY :: ========>",o_tot_ene 00794 ! Debug Potential 00795 dq = 0.001_dp 00796 tot_ene = 0.0_dp 00797 DO i = 1, nparticles 00798 DO k = 1, 2 00799 lcharges = charges 00800 lcharges(i) = charges(i) + (-1.0_dp)**k * dq 00801 forces_glob = 0.0_dp 00802 forces_local= 0.0_dp 00803 pv_local = 0.0_dp 00804 pv_glob = 0.0_dp 00805 energy_glob = 0.0_dp 00806 energy_local= 0.0_dp 00807 e_neut = 0.0_dp 00808 e_self = 0.0_dp 00809 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set,& 00810 local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .FALSE., .FALSE.,& 00811 .FALSE., radii, lcharges, dipoles, quadrupoles, iw=iw, do_debug=.FALSE., error=error ) 00812 ene(k) = energy_local + energy_glob + e_neut + e_self 00813 END DO 00814 pot = (ene(2)-ene(1))/(2.0_dp*dq) 00815 WRITE(*,'(A,I8,3(A,F15.9))')"POTENTIAL FOR ATOM: ",i, " NUMERICAL: ",pot, " ANALYTICAL: ",efield0(i),& 00816 " ERROR: ",pot-efield0(i) 00817 tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i) 00818 END DO 00819 WRITE(*,*)"ENERGIES: ",o_tot_ene, tot_ene, o_tot_ene-tot_ene 00820 WRITE(*,'(/,/,/)') 00821 ! Debug Field 00822 dq = 0.001_dp 00823 DO i = 1, nparticles 00824 coord = particle_set(i)%r 00825 DO j = 1, 3 00826 DO k = 1, 2 00827 particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k * dq 00828 00829 ! Rebuild neighbor lists 00830 CALL list_control ( atomic_kind_set, particle_set, local_particles, & 00831 cell, nonbond_env, logger%para_env, mm_section, & 00832 shell_particle_set, core_particle_set, error=error) 00833 00834 forces_glob = 0.0_dp 00835 forces_local= 0.0_dp 00836 pv_local = 0.0_dp 00837 pv_glob = 0.0_dp 00838 energy_glob = 0.0_dp 00839 energy_local= 0.0_dp 00840 e_neut = 0.0_dp 00841 e_self = 0.0_dp 00842 efield0 = 0.0_dp 00843 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set,& 00844 local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .TRUE., .TRUE.,& 00845 .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob,& 00846 efield0, iw=iw, do_debug=.FALSE., error=error) 00847 ene(k) = efield0(i) 00848 particle_set(i)%r(j) = coord(j) 00849 END DO 00850 efield1n(j) = -(ene(2)-ene(1))/(2.0_dp*dq) 00851 END DO 00852 WRITE(*,'(/,A,I8)')"FIELD FOR ATOM: ",i 00853 WRITE(*,'(A,3F15.9)')" NUMERICAL: ",efield1n, " ANALYTICAL: ",efield1(:,i),& 00854 " ERROR: ",efield1n-efield1(:,i) 00855 IF (task(2)) THEN 00856 tot_ene = tot_ene - 0.5_dp*DOT_PRODUCT(efield1(:,i),dipoles(:,i)) 00857 END IF 00858 END DO 00859 WRITE(*,*)"ENERGIES: ",o_tot_ene, tot_ene, o_tot_ene-tot_ene 00860 00861 ! Debug Field Gradient 00862 dq = 0.0001_dp 00863 DO i = 1, nparticles 00864 coord = particle_set(i)%r 00865 DO j = 1, 3 00866 DO k = 1, 2 00867 particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k * dq 00868 00869 ! Rebuild neighbor lists 00870 CALL list_control ( atomic_kind_set, particle_set, local_particles, & 00871 cell, nonbond_env, logger%para_env, mm_section, & 00872 shell_particle_set, core_particle_set, error=error) 00873 00874 forces_glob = 0.0_dp 00875 forces_local= 0.0_dp 00876 pv_local = 0.0_dp 00877 pv_glob = 0.0_dp 00878 energy_glob = 0.0_dp 00879 energy_local= 0.0_dp 00880 e_neut = 0.0_dp 00881 e_self = 0.0_dp 00882 efield1 = 0.0_dp 00883 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set,& 00884 local_particles, energy_local, energy_glob, e_neut, e_self, task,.FALSE., .TRUE., .TRUE.,& 00885 .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob,& 00886 efield1=efield1, iw=iw, do_debug=.FALSE., error=error) 00887 enev(:,k) = efield1(:,i) 00888 particle_set(i)%r(j) = coord(j) 00889 END DO 00890 efield2n(:,j) = (enev(:,2)-enev(:,1))/(2.0_dp*dq) 00891 END DO 00892 WRITE(*,'(/,A,I8)')"FIELD GRADIENT FOR ATOM: ",i 00893 WRITE(*,'(A,9F15.9)')" NUMERICAL: ",efield2n,& 00894 " ANALYTICAL: ",efield2(:,i),& 00895 " ERROR: ",RESHAPE(efield2n,(/9/))-efield2(:,i) 00896 END DO 00897 END SUBROUTINE debug_ewald_multipoles_fields 00898 00899 ! ***************************************************************************** 00904 SUBROUTINE debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, cell,& 00905 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw, & 00906 error ) 00907 USE cell_types, ONLY: cell_type 00908 USE distribution_1d_types, ONLY: distribution_1d_type 00909 USE ewald_environment_types, ONLY: ewald_environment_type 00910 USE ewald_pw_types, ONLY: ewald_pw_type 00911 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 00912 neighbor_kind_pairs_type 00913 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 00914 fist_nonbond_env_type,& 00915 pos_type 00916 USE kinds, ONLY: dp 00917 USE particle_types, ONLY: particle_type 00918 USE ewalds_multipole, ONLY: ewald_multipole_evaluate 00919 USE fist_neighbor_list_control, ONLY: list_control 00920 #include "cp_common_uses.h" 00921 IMPLICIT NONE 00922 TYPE(ewald_environment_type), POINTER :: ewald_env 00923 TYPE(ewald_pw_type), POINTER :: ewald_pw 00924 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 00925 TYPE(cell_type), POINTER :: cell 00926 TYPE(particle_type), POINTER :: particle_set(:) 00927 TYPE(distribution_1d_type), POINTER :: local_particles 00928 REAL(KIND=dp), DIMENSION(:), 00929 POINTER, OPTIONAL :: radii, charges 00930 REAL(KIND=dp), DIMENSION(:, :), 00931 POINTER, OPTIONAL :: dipoles 00932 REAL(KIND=dp), DIMENSION(:, :, :), 00933 POINTER, OPTIONAL :: quadrupoles 00934 LOGICAL, DIMENSION(3), INTENT(IN) :: task 00935 INTEGER, INTENT(IN) :: iw 00936 TYPE(cp_error_type), INTENT(inout) :: error 00937 00938 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipoles_fields', 00939 routineP = "ewalds_multipole_debug"//':'//routineN 00940 00941 REAL(KIND=dp), POINTER, DIMENSION(:) :: efield0 00942 REAL(KIND=dp), POINTER, DIMENSION(:,:) :: ldipoles 00943 REAL(KIND=dp), POINTER, DIMENSION(:,:) :: lquadrupoles 00944 REAL(KIND=dp) :: dq, dr, energy_local, energy_glob, e_neut, 00945 e_self, pv_local(3,3), pv_glob(3,3), ene(2), pot, coord(3), efield1n(3), tot_ene, 00946 o_tot_ene, enev(3,2), efield2n(3,3), prod 00947 REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE 00948 :: forces_local, forces_glob, efield1, efield2 00949 INTEGER :: nparticles,i,k,j,nparticle_local,iparticle_kind,ind 00950 TYPE(particle_type), POINTER, DIMENSION(:) :: shell_particle_set, core_particle_set 00951 TYPE(cp_logger_type), POINTER :: logger 00952 00953 NULLIFY(ldipoles, lquadrupoles, shell_particle_set, core_particle_set) 00954 NULLIFY(logger) 00955 logger => cp_error_get_logger(error) 00956 00957 nparticles = SIZE(particle_set) 00958 nparticle_local = 0 00959 DO iparticle_kind=1,SIZE(local_particles%n_el) 00960 nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind) 00961 END DO 00962 ALLOCATE(forces_glob(3,nparticles)) 00963 ALLOCATE(forces_local(3,nparticle_local)) 00964 ALLOCATE(efield0(nparticles)) 00965 ALLOCATE(efield1(3,nparticles)) 00966 ALLOCATE(efield2(9,nparticles)) 00967 forces_glob = 0.0_dp 00968 forces_local= 0.0_dp 00969 efield0 = 0.0_dp 00970 efield1 = 0.0_dp 00971 efield2 = 0.0_dp 00972 pv_local = 0.0_dp 00973 pv_glob = 0.0_dp 00974 energy_glob = 0.0_dp 00975 energy_local= 0.0_dp 00976 e_neut = 0.0_dp 00977 e_self = 0.0_dp 00978 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set,& 00979 local_particles, energy_local, energy_glob, e_neut, e_self, task,.FALSE.,.TRUE., .TRUE.,& 00980 .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob,& 00981 efield0, efield1, efield2, iw, do_debug=.FALSE., error=error ) 00982 o_tot_ene = energy_local + energy_glob + e_neut + e_self 00983 WRITE(*,*)"TOTAL ENERGY :: ========>",o_tot_ene 00984 00985 ! Debug Potential 00986 tot_ene = 0.0_dp 00987 IF (task(1)) THEN 00988 DO i = 1, nparticles 00989 tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i) 00990 END DO 00991 WRITE(*,*)"ENERGIES: ",o_tot_ene, tot_ene, o_tot_ene-tot_ene 00992 WRITE(*,'(/,/,/)') 00993 END IF 00994 00995 ! Debug Field 00996 IF (task(2)) THEN 00997 DO i = 1, nparticles 00998 tot_ene = tot_ene - 0.5_dp*DOT_PRODUCT(efield1(:,i),dipoles(:,i)) 00999 END DO 01000 WRITE(*,*)"ENERGIES: ",o_tot_ene, tot_ene, o_tot_ene-tot_ene 01001 WRITE(*,'(/,/,/)') 01002 END IF 01003 01004 ! Debug Field Gradient 01005 IF (task(3)) THEN 01006 DO i = 1, nparticles 01007 ind = 0 01008 prod = 0.0_dp 01009 DO j = 1,3 01010 DO k = 1,3 01011 ind = ind + 1 01012 prod = prod + efield2(ind,i)*quadrupoles(j,k,i) 01013 END DO 01014 END DO 01015 tot_ene = tot_ene - 0.5_dp*(1.0_dp/3.0_dp)*prod 01016 END DO 01017 WRITE(*,*)"ENERGIES: ",o_tot_ene, tot_ene, o_tot_ene-tot_ene 01018 WRITE(*,'(/,/,/)') 01019 END IF 01020 01021 END SUBROUTINE debug_ewald_multipoles_fields2
1.7.3