|
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 ! ***************************************************************************** 00015 MODULE qs_linres_current 00016 00017 USE array_types, ONLY: array_i1d_obj,& 00018 array_new,& 00019 array_nullify,& 00020 array_release 00021 USE atomic_kind_types, ONLY: atomic_kind_type,& 00022 get_atomic_kind,& 00023 get_atomic_kind_set 00024 USE basis_set_types, ONLY: get_gto_basis_set,& 00025 gto_basis_set_p_type,& 00026 gto_basis_set_type 00027 USE cell_types, ONLY: cell_type,& 00028 pbc 00029 USE cp_array_i_utils, ONLY: cp_2d_i_p_type 00030 USE cp_array_r_utils, ONLY: cp_2d_r_p_type 00031 USE cp_control_types, ONLY: dft_control_type 00032 USE cp_dbcsr_interface, ONLY: & 00033 cp_dbcsr_col_block_sizes, cp_dbcsr_copy, cp_dbcsr_create, & 00034 cp_dbcsr_distribution, cp_dbcsr_finalize, cp_dbcsr_get_block_p, & 00035 cp_dbcsr_get_data_size, cp_dbcsr_get_matrix_type, & 00036 cp_dbcsr_get_num_blocks, cp_dbcsr_init, cp_dbcsr_row_block_sizes, & 00037 cp_dbcsr_set 00038 USE cp_dbcsr_operations, ONLY: cp_dbcsr_add_block_node,& 00039 cp_dbcsr_alloc_block_from_nbl,& 00040 cp_dbcsr_allocate_matrix_set,& 00041 cp_dbcsr_deallocate_matrix,& 00042 cp_dbcsr_deallocate_matrix_set,& 00043 cp_dbcsr_plus_fm_fm_t,& 00044 cp_dbcsr_sm_fm_multiply 00045 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00046 cp_dbcsr_type 00047 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& 00048 cp_fm_trace 00049 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00050 cp_fm_struct_release,& 00051 cp_fm_struct_type 00052 USE cp_fm_types, ONLY: cp_fm_create,& 00053 cp_fm_p_type,& 00054 cp_fm_release,& 00055 cp_fm_set_all,& 00056 cp_fm_to_fm,& 00057 cp_fm_type 00058 USE cp_output_handling, ONLY: cp_p_file,& 00059 cp_print_key_finished_output,& 00060 cp_print_key_should_output,& 00061 cp_print_key_unit_nr 00062 USE cp_para_types, ONLY: cp_para_env_type 00063 USE cp_subsys_types, ONLY: cp_subsys_get,& 00064 cp_subsys_type 00065 USE cube_utils, ONLY: cube_info_type 00066 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00067 dbcsr_type_antisymmetric,& 00068 dbcsr_type_no_symmetry 00069 USE dbcsr_util, ONLY: convert_offsets_to_sizes 00070 USE gaussian_gridlevels, ONLY: gridlevel_info_type 00071 USE input_constants, ONLY: current_gauge_atom 00072 USE input_section_types, ONLY: section_get_ivals,& 00073 section_get_lval,& 00074 section_vals_get_subs_vals,& 00075 section_vals_type 00076 USE kinds, ONLY: default_path_length,& 00077 dp,& 00078 int_8 00079 USE lgrid_types, ONLY: lgrid_type 00080 USE mathconstants, ONLY: twopi 00081 USE memory_utilities, ONLY: reallocate 00082 USE message_passing 00083 USE orbital_pointers, ONLY: ncoset 00084 USE particle_list_types, ONLY: particle_list_type 00085 USE particle_types, ONLY: get_particle_set,& 00086 particle_type 00087 USE pw_env_types, ONLY: pw_env_get,& 00088 pw_env_type 00089 USE pw_methods, ONLY: pw_axpy,& 00090 pw_integrate_function,& 00091 pw_scale,& 00092 pw_zero 00093 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00094 pw_pool_give_back_pw,& 00095 pw_pool_type 00096 USE pw_types, ONLY: REALDATA3D,& 00097 REALSPACE,& 00098 pw_p_type 00099 USE qs_collocate_density, ONLY: collocate_pgf_product_rspace,& 00100 density_rs2pw 00101 USE qs_environment_types, ONLY: get_qs_env,& 00102 qs_environment_type 00103 USE qs_linres_atom_current, ONLY: calculate_jrho_atom,& 00104 calculate_jrho_atom_coeff,& 00105 calculate_jrho_atom_rad 00106 USE qs_linres_op, ONLY: fac_vecp,& 00107 fm_scale_by_pbc_AC,& 00108 ind_m2,& 00109 set_vecp,& 00110 set_vecp_rev 00111 USE qs_linres_types, ONLY: current_env_type,& 00112 get_current_env,& 00113 realspaces_grid_p_type 00114 USE qs_matrix_pools, ONLY: qs_matrix_pools_type 00115 USE qs_mo_types, ONLY: get_mo_set,& 00116 mo_set_p_type 00117 USE qs_modify_pab_block, ONLY: FUNC_AB,& 00118 FUNC_ADBmDAB,& 00119 FUNC_ARDBmDARB 00120 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 00121 neighbor_list_iterate,& 00122 neighbor_list_iterator_create,& 00123 neighbor_list_iterator_p_type,& 00124 neighbor_list_iterator_release,& 00125 neighbor_list_set_p_type 00126 USE qs_operators_ao, ONLY: build_lin_mom_matrix,& 00127 rRc_xyz_der_ao 00128 USE realspace_grid_cube, ONLY: pw_to_cube 00129 USE realspace_grid_types, ONLY: & 00130 realspace_grid_desc_p_type, realspace_grid_desc_type, & 00131 realspace_grid_p_type, realspace_grid_type, rs_grid_create, & 00132 rs_grid_mult_and_add, rs_grid_release, rs_grid_retain, rs_grid_zero 00133 USE task_list_methods, ONLY: distribute_tasks,& 00134 int2pair,& 00135 rs_distribute_matrix,& 00136 task_list_inner_loop 00137 USE termination, ONLY: stop_program 00138 USE timings, ONLY: timeset,& 00139 timestop 00140 #include "cp_common_uses.h" 00141 00142 IMPLICIT NONE 00143 00144 PRIVATE 00145 00146 ! *** Public subroutines *** 00147 PUBLIC :: current_build_current,calculate_jrho_resp,current_build_chi,current_set_gauge 00148 00149 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_current' 00150 00151 TYPE box_type 00152 INTEGER :: n 00153 REAL(dp), POINTER, DIMENSION(:,:) :: r 00154 END TYPE box_type 00155 REAL(dp), DIMENSION(3, 3, 3),PARAMETER :: Levi_Civita = RESHAPE((/ 00156 0.0_dp, 0.0_dp,0.0_dp,0.0_dp,0.0_dp,-1.0_dp, 0.0_dp,1.0_dp,0.0_dp, 00157 0.0_dp, 0.0_dp,1.0_dp,0.0_dp,0.0_dp, 0.0_dp,-1.0_dp,0.0_dp,0.0_dp, 00158 0.0_dp,-1.0_dp,0.0_dp,1.0_dp,0.0_dp, 0.0_dp, 0.0_dp,0.0_dp,0.0_dp/),(/3,3,3/)) 00159 00160 00161 CONTAINS 00162 00163 ! ***************************************************************************** 00186 SUBROUTINE current_build_current(current_env,qs_env,iB,error) 00187 ! 00188 TYPE(current_env_type) :: current_env 00189 TYPE(qs_environment_type), POINTER :: qs_env 00190 INTEGER, INTENT(IN) :: iB 00191 TYPE(cp_error_type), INTENT(INOUT) :: error 00192 00193 CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_current', 00194 routineP = moduleN//':'//routineN 00195 00196 CHARACTER(LEN=default_path_length) :: ext, filename, my_pos 00197 INTEGER :: handle, idir, iiB, iiiB, ispin, istat, istate, j, jstate, nao, 00198 natom, nmo, nspins, nstates(2), output_unit, unit_nr 00199 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 00200 INTEGER, DIMENSION(:), POINTER :: rbs 00201 LOGICAL :: append_cube, failure, gapw 00202 REAL(dp) :: dk(3), jrho_tot_G(3,3), 00203 jrho_tot_R(3,3), maxocc, 00204 scale_fac 00205 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ddk 00206 REAL(dp), EXTERNAL :: DDOT 00207 TYPE(array_i1d_obj) :: row_blk_sizes 00208 TYPE(cell_type), POINTER :: cell 00209 TYPE(cp_2d_i_p_type), DIMENSION(:), 00210 POINTER :: center_list 00211 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00212 POINTER :: density_matrix0, 00213 density_matrix_a, 00214 density_matrix_ii, 00215 density_matrix_iii 00216 TYPE(cp_fm_p_type), DIMENSION(:), 00217 POINTER :: p_psi1, psi0_order, psi1 00218 TYPE(cp_fm_p_type), DIMENSION(:, :), 00219 POINTER :: psi1_D, psi1_p, psi1_rxp 00220 TYPE(cp_fm_type), POINTER :: mo_coeff, psi_a_iB, psi_buf 00221 TYPE(cp_logger_type), POINTER :: logger 00222 TYPE(cp_para_env_type), POINTER :: para_env 00223 TYPE(cp_subsys_type), POINTER :: subsys 00224 TYPE(dbcsr_distribution_obj), POINTER :: dbcsr_dist 00225 TYPE(dft_control_type), POINTER :: dft_control 00226 TYPE(mo_set_p_type), DIMENSION(:), 00227 POINTER :: mos 00228 TYPE(neighbor_list_set_p_type), 00229 DIMENSION(:), POINTER :: sab_all 00230 TYPE(particle_list_type), POINTER :: particles 00231 TYPE(particle_type), DIMENSION(:), 00232 POINTER :: particle_set 00233 TYPE(pw_env_type), POINTER :: pw_env 00234 TYPE(pw_p_type) :: wf_r 00235 TYPE(pw_p_type), POINTER :: jrho_gspace, jrho_rspace 00236 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00237 TYPE(qs_matrix_pools_type), POINTER :: mpools 00238 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00239 TYPE(section_vals_type), POINTER :: current_section 00240 00241 failure = .FALSE. 00242 CALL timeset(routineN,handle) 00243 ! 00244 NULLIFY(jrho_rspace,jrho_gspace,logger,current_section,density_matrix0,density_matrix_a,& 00245 & density_matrix_ii,density_matrix_iii,cell,dft_control,mos,& 00246 & particle_set,pw_env,auxbas_rs_desc,auxbas_pw_pool,& 00247 & para_env,center_list,mo_coeff,psi_a_iB,& 00248 & psi1,p_psi1,psi1_p,psi1_D,psi1_rxp,psi0_order,sab_all) 00249 00250 logger => cp_error_get_logger(error) 00251 output_unit= cp_logger_get_default_io_unit(logger) 00252 ! 00253 ! 00254 CALL get_current_env(current_env=current_env,& 00255 & center_list=center_list,& 00256 & psi1_rxp=psi1_rxp,& 00257 & psi1_D=psi1_D,& 00258 & psi1_p=psi1_p,& 00259 & psi0_order=psi0_order,& 00260 & nstates=nstates,& 00261 & nao=nao,& 00262 & error=error) 00263 ! 00264 ! 00265 CALL get_qs_env(qs_env=qs_env,& 00266 & cell=cell,& 00267 & dft_control=dft_control,& 00268 & mos=mos,& 00269 & mpools=mpools,& 00270 & pw_env=pw_env,& 00271 & para_env=para_env,& 00272 & subsys=subsys,& 00273 & sab_all=sab_all,& 00274 & particle_set=particle_set,& 00275 & dbcsr_dist=dbcsr_dist,& 00276 & error=error) 00277 00278 CALL cp_subsys_get(subsys,particles=particles,error=error) 00279 00280 gapw = dft_control%qs_control%gapw 00281 nspins = dft_control%nspins 00282 natom = SIZE(particle_set,1) 00283 ! 00284 ! allocate temporary arrays 00285 ALLOCATE(psi1(nspins),p_psi1(nspins),STAT=istat) 00286 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00287 DO ispin = 1,nspins 00288 CALL cp_fm_create( psi1(ispin)%matrix,psi0_order(ispin)%matrix%matrix_struct,error=error) 00289 CALL cp_fm_create(p_psi1(ispin)%matrix,psi0_order(ispin)%matrix%matrix_struct,error=error) 00290 CALL cp_fm_set_all( psi1(ispin)%matrix,0.0_dp,error=error) 00291 CALL cp_fm_set_all(p_psi1(ispin)%matrix,0.0_dp,error=error) 00292 ENDDO 00293 ! 00294 ! 00295 CALL cp_dbcsr_allocate_matrix_set(density_matrix0,nspins,error=error) 00296 CALL cp_dbcsr_allocate_matrix_set(density_matrix_a,nspins,error=error) 00297 CALL cp_dbcsr_allocate_matrix_set(density_matrix_ii,nspins,error=error) 00298 CALL cp_dbcsr_allocate_matrix_set(density_matrix_iii,nspins,error=error) 00299 ! 00300 ! prepare for allocation 00301 ALLOCATE (first_sgf(natom),STAT=istat) 00302 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00303 ALLOCATE (last_sgf(natom),STAT=istat) 00304 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00305 CALL get_particle_set(particle_set=particle_set,& 00306 first_sgf=first_sgf,& 00307 last_sgf=last_sgf,error=error) 00308 ALLOCATE (rbs(natom), STAT=istat) 00309 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00310 CALL convert_offsets_to_sizes (first_sgf, rbs, last_sgf) 00311 CALL array_nullify (row_blk_sizes) 00312 CALL array_new (row_blk_sizes, rbs, gift=.TRUE.) 00313 DEALLOCATE (first_sgf,STAT=istat) 00314 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00315 DEALLOCATE (last_sgf,STAT=istat) 00316 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00317 ! 00318 ! 00319 DO ispin = 1,nspins 00320 ! 00321 !density_matrix0 00322 ALLOCATE(density_matrix0(ispin)%matrix) 00323 CALL cp_dbcsr_init(density_matrix0(ispin)%matrix,error=error) 00324 CALL cp_dbcsr_create(matrix=density_matrix0(ispin)%matrix, & 00325 name="density_matrix0", & 00326 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry,& 00327 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00328 nblks=0, nze=0, mutable_work=.TRUE., & 00329 error=error) 00330 CALL cp_dbcsr_alloc_block_from_nbl(density_matrix0(ispin)%matrix,sab_all,error=error) 00331 ! 00332 !density_matrix_a 00333 ALLOCATE(density_matrix_a(ispin)%matrix) 00334 CALL cp_dbcsr_init(density_matrix_a(ispin)%matrix,error=error) 00335 CALL cp_dbcsr_copy(density_matrix_a(ispin)%matrix,density_matrix0(ispin)%matrix,& 00336 name="density_matrix_a",error=error) 00337 ! 00338 !density_matrix_ii 00339 ALLOCATE(density_matrix_ii(ispin)%matrix) 00340 CALL cp_dbcsr_init(density_matrix_ii(ispin)%matrix,error=error) 00341 CALL cp_dbcsr_copy(density_matrix_ii(ispin)%matrix,density_matrix0(ispin)%matrix,& 00342 name="density_matrix_ii",error=error) 00343 ! 00344 !density_matrix_iii 00345 ALLOCATE(density_matrix_iii(ispin)%matrix) 00346 CALL cp_dbcsr_init(density_matrix_iii(ispin)%matrix,error=error) 00347 CALL cp_dbcsr_copy(density_matrix_iii(ispin)%matrix,density_matrix0(ispin)%matrix,& 00348 name="density_matrix_iii",error=error) 00349 ENDDO 00350 ! 00351 CALL array_release (row_blk_sizes) 00352 ! 00353 ! 00354 current_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES%CURRENT",error=error) 00355 ! 00356 IF(.NOT. failure) THEN 00357 ! 00358 ! 00359 jrho_tot_G = 0.0_dp 00360 jrho_tot_R = 0.0_dp 00361 ! 00362 ! Lets go! 00363 CALL set_vecp(iB,iiB,iiiB) 00364 DO ispin = 1,nspins 00365 nmo = nstates(ispin) 00366 mo_coeff => psi0_order(ispin)%matrix 00367 !maxocc = max_occ(ispin) 00368 ! 00369 CALL get_mo_set(mo_set=mos(ispin)%mo_set,maxocc=maxocc) 00370 ! 00371 ! 00372 ! Build the first density matrix 00373 CALL cp_dbcsr_set(density_matrix0(ispin)%matrix,0.0_dp,error=error) 00374 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix0(ispin)%matrix ,& 00375 & matrix_v=mo_coeff,matrix_g=mo_coeff,& 00376 & ncol=nmo,alpha=maxocc,error=error) 00377 ! 00378 ! Allocate buffer vectors 00379 ALLOCATE(ddk(3,nmo),STAT=istat) 00380 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00381 ! 00382 ! Construct the 3 density matrices for the field in direction iB 00383 ! 00384 ! First the full matrix psi_a_iB 00385 psi_a_iB => psi1(ispin)%matrix 00386 psi_buf => p_psi1(ispin)%matrix 00387 CALL cp_fm_set_all(psi_a_iB,0.0_dp,error=error) 00388 CALL cp_fm_set_all(psi_buf,0.0_dp,error=error) 00389 ! psi_a_iB = - (R_\nu-dk)_ii psi1_piiiB + (R_\nu-dk)_iii psi1_piiB 00390 ! 00391 ! contributions from the response psi1_p_ii and psi1_p_iii 00392 DO istate = 1,current_env%nbr_center(ispin) 00393 dk(1:3) = current_env%centers_set(ispin)%array(1:3,istate) 00394 ! 00395 ! Copy the vector in the full matrix psi1 00396 !nstate_loc = center_list(ispin)%array(1,icenter+1)-center_list(ispin)%array(1,icenter) 00397 DO j = center_list(ispin)%array(1,istate),center_list(ispin)%array(1,istate+1)-1 00398 jstate = center_list(ispin)%array(2,j) 00399 CALL cp_fm_to_fm(psi1_p(ispin,iiB)%matrix,psi_a_iB,1,jstate,jstate) 00400 CALL cp_fm_to_fm(psi1_p(ispin,iiiB)%matrix,psi_buf,1,jstate,jstate) 00401 ddk(:,jstate) = dk(1:3) 00402 ENDDO 00403 ENDDO ! istate 00404 CALL fm_scale_by_pbc_AC(psi_a_iB,current_env%basisfun_center,ddk,cell,iiiB,error) 00405 CALL fm_scale_by_pbc_AC(psi_buf,current_env%basisfun_center,ddk,cell,iiB,error) 00406 CALL cp_fm_scale_and_add(-1.0_dp,psi_a_iB,1.0_dp,psi_buf,error=error) 00407 ! 00408 !psi_a_iB = psi_a_iB + psi1_rxp 00409 ! 00410 ! contribution from the response psi1_rxp 00411 CALL cp_fm_scale_and_add(-1.0_dp,psi_a_iB,1.0_dp,psi1_rxp(ispin,iB)%matrix,error=error) 00412 ! 00413 !psi_a_iB = psi_a_iB - psi1_D 00414 IF(current_env%full) THEN 00415 ! 00416 ! contribution from the response psi1_D 00417 CALL cp_fm_scale_and_add(1.0_dp,psi_a_iB,-1.0_dp,psi1_D(ispin,iB)%matrix,error=error) 00418 ENDIF 00419 ! 00420 ! Multiply by the occupation number for the density matrix 00421 ! 00422 ! Build the first density matrix 00423 CALL cp_dbcsr_set(density_matrix_a(ispin)%matrix ,0.0_dp,error=error) 00424 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_a(ispin)%matrix ,& 00425 & matrix_v=mo_coeff,matrix_g=psi_a_iB,& 00426 & ncol=nmo,alpha=maxocc,error=error) 00427 ! 00428 ! Build the second density matrix 00429 CALL cp_dbcsr_set(density_matrix_iii(ispin)%matrix,0.0_dp,error=error) 00430 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_iii(ispin)%matrix ,& 00431 & matrix_v=mo_coeff,matrix_g=psi1_p(ispin,iiiB)%matrix,& 00432 & ncol=nmo,alpha=maxocc,error=error) 00433 ! 00434 ! Build the third density matrix 00435 CALL cp_dbcsr_set(density_matrix_ii(ispin)%matrix,0.0_dp,error=error) 00436 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_ii(ispin)%matrix ,& 00437 & matrix_v=mo_coeff,matrix_g=psi1_p(ispin,iiB)%matrix,& 00438 & ncol=nmo,alpha=maxocc,error=error) 00439 DO idir = 1,3 00440 ! 00441 ! Calculate the current density on the pw grid (only soft if GAPW) 00442 ! idir is the cartesian component of the response current density 00443 ! generated by the magnetic field pointing in cartesian direction iB 00444 ! Use the qs_rho_type already used for rho during the scf 00445 jrho_rspace => current_env%jrho1_set(idir)%rho%rho_r(ispin) 00446 jrho_gspace => current_env%jrho1_set(idir)%rho%rho_g(ispin) 00447 CALL pw_zero(jrho_rspace%pw,error=error) 00448 CALL pw_zero(jrho_gspace%pw,error=error) 00449 CALL calculate_jrho_resp(density_matrix0(ispin)%matrix, & 00450 & density_matrix_a(ispin)%matrix, & 00451 & density_matrix_ii(ispin)%matrix, & 00452 & density_matrix_iii(ispin)%matrix, & 00453 & iB,idir,jrho_rspace,jrho_gspace,qs_env, & 00454 & current_env,gapw,error=error) 00455 00456 scale_fac = cell%deth / twopi 00457 CALL pw_scale(jrho_rspace%pw,scale_fac,error=error) 00458 CALL pw_scale(jrho_gspace%pw,scale_fac,error=error) 00459 00460 jrho_tot_G(idir,iB) = pw_integrate_function(jrho_gspace%pw,isign=-1,error=error) 00461 jrho_tot_R(idir,iB) = pw_integrate_function(jrho_rspace%pw,isign=-1,error=error) 00462 00463 IF(output_unit>0) THEN 00464 WRITE(output_unit,'(T2,2(A,E24.16))') 'Integrated j_'& 00465 &//ACHAR(idir+119)//ACHAR(iB+119)//'(r): G-space=',& 00466 &jrho_tot_G(idir,iB),' R-space=',jrho_tot_R(idir,iB) 00467 ENDIF 00468 ! 00469 ENDDO ! idir 00470 ! 00471 ! Deallocate buffer vectors 00472 DEALLOCATE(ddk,STAT=istat) 00473 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00474 ! 00475 ENDDO ! ispin 00476 00477 IF(gapw) THEN 00478 DO idir = 1,3 00479 ! 00480 ! compute the atomic response current densities on the spherical grids 00481 ! First the sparse matrices are multiplied by the expansion coefficients 00482 ! this is the equivalent of the CPC for the charge density 00483 CALL calculate_jrho_atom_coeff(qs_env,current_env, & 00484 & density_matrix0, & 00485 & density_matrix_a, & 00486 & density_matrix_ii, & 00487 & density_matrix_iii, & 00488 & iB,idir,error=error) 00489 ! 00490 ! Then the radial parts are computed on the local radial grid, atom by atom 00491 ! 8 functions are computed for each atom, per grid point 00492 ! and per LM angular momentum. The multiplication by the Clebsh-Gordon 00493 ! coefficients or they correspondent for the derivatives, is also done here 00494 CALL calculate_jrho_atom_rad(qs_env,current_env,idir,error=error) 00495 ! 00496 ! The current on the atomic grids 00497 CALL calculate_jrho_atom(current_env,qs_env,iB,idir,error=error) 00498 ENDDO ! idir 00499 ENDIF 00500 ! 00501 ! Cube files 00502 IF(BTEST(cp_print_key_should_output(logger%iter_info,current_section,& 00503 & "PRINT%CURRENT_CUBES",error=error),cp_p_file)) THEN 00504 append_cube = section_get_lval(current_section,"PRINT%CURRENT_CUBES%APPEND",error=error) 00505 my_pos = "REWIND" 00506 IF(append_cube) THEN 00507 my_pos = "APPEND" 00508 END IF 00509 ! 00510 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc,& 00511 & auxbas_pw_pool=auxbas_pw_pool,& 00512 & error=error) 00513 ! 00514 CALL pw_pool_create_pw(auxbas_pw_pool,wf_r%pw,use_data=REALDATA3D,& 00515 & in_space=REALSPACE,error=error) 00516 ! 00517 DO idir = 1,3 00518 CALL pw_zero(wf_r%pw,error=error) 00519 ! 00520 DO ispin =1 ,nspins 00521 CALL pw_axpy(current_env%jrho1_set(idir)%rho%rho_r(ispin)%pw,wf_r%pw,1.0_dp,& 00522 error=error) 00523 ENDDO 00524 ! 00525 IF(gapw) THEN 00526 ! Add the local hard and soft contributions 00527 ! This can be done atom by atom by a spline extrapolation of the values 00528 ! on the spherical grid to the grid points. 00529 CALL stop_program(routineN,moduleN,__LINE__,"GAPW needs to be finalized") 00530 ENDIF 00531 filename="jresp" 00532 WRITE(ext,'(a2,I1,a2,I1,a5)') "iB",iB,"_d",idir,".cube" 00533 WRITE(ext,'(a2,a1,a2,a1,a5)') "iB",ACHAR(iB+119),"_d",ACHAR(idir+119),".cube" 00534 unit_nr=cp_print_key_unit_nr(logger,current_section,"PRINT%CURRENT_CUBES",& 00535 & extension=TRIM(ext),middle_name=TRIM(filename),& 00536 & log_filename=.FALSE.,file_position=my_pos,& 00537 & error=error) 00538 00539 CALL pw_to_cube(wf_r%pw,unit_nr,"RESPONSE CURRENT DENSITY ",& 00540 & particles=particles,& 00541 & stride=section_get_ivals(current_section,"PRINT%CURRENT_CUBES%STRIDE",& 00542 & error=error),error=error) 00543 CALL cp_print_key_finished_output(unit_nr,logger,current_section,& 00544 & "PRINT%CURRENT_CUBES",error=error) 00545 ENDDO 00546 ! 00547 CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_r%pw,error=error) 00548 ENDIF ! current cube 00549 ! 00550 ! Integrated current response checksum 00551 IF(output_unit>0) THEN 00552 WRITE(output_unit,'(T2,A,E24.16)') 'CheckSum R-integrated j=',& 00553 & SQRT(DDOT(9,jrho_tot_R(1,1),1,jrho_tot_R(1,1),1)) 00554 ENDIF 00555 ! 00556 ! 00557 ! Dellocate grids for the calculation of jrho and the shift 00558 DO ispin = 1,nspins 00559 CALL cp_fm_release( psi1(ispin)%matrix,error) 00560 CALL cp_fm_release(p_psi1(ispin)%matrix,error) 00561 ENDDO 00562 DEALLOCATE(psi1,p_psi1,STAT=istat) 00563 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00564 00565 CALL cp_dbcsr_deallocate_matrix_set(density_matrix0,error=error) 00566 CALL cp_dbcsr_deallocate_matrix_set(density_matrix_a,error=error) 00567 CALL cp_dbcsr_deallocate_matrix_set(density_matrix_ii,error=error) 00568 CALL cp_dbcsr_deallocate_matrix_set(density_matrix_iii,error=error) 00569 00570 ! 00571 ENDIF ! failure 00572 ! 00573 ! Finalize 00574 CALL timestop(handle) 00575 ! 00576 END SUBROUTINE current_build_current 00577 00578 ! ***************************************************************************** 00605 SUBROUTINE calculate_jrho_resp(mat_d0,mat_jp,mat_jp_rii,mat_jp_riii,iB,idir,& 00606 current_rs, current_gs, qs_env, current_env, soft_valid, error) 00607 00608 TYPE(cp_dbcsr_type), POINTER :: mat_d0, mat_jp, mat_jp_rii, 00609 mat_jp_riii 00610 INTEGER, INTENT(IN) :: iB, idir 00611 TYPE(pw_p_type), INTENT(INOUT) :: current_rs, current_gs 00612 TYPE(qs_environment_type), POINTER :: qs_env 00613 TYPE(current_env_type) :: current_env 00614 LOGICAL, INTENT(IN), OPTIONAL :: soft_valid 00615 TYPE(cp_error_type), INTENT(inout) :: error 00616 00617 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_jrho_resp', 00618 routineP = moduleN//':'//routineN 00619 INTEGER, PARAMETER :: add_tasks = 1000, 00620 max_tasks = 2000 00621 REAL(kind=dp), PARAMETER :: mult_tasks = 2.0_dp 00622 00623 INTEGER :: bcol, brow, curr_tasks, handle, i, iatom, iatom_old, idir2, 00624 igrid_level, iiB, iiiB, ikind, ikind_old, ipgf, iset, iset_old, istat, 00625 itask, ithread, jatom, jatom_old, jkind, jkind_old, jpgf, jset, 00626 jset_old, maxco, maxpgf, maxset, maxsgf, maxsgf_set, na1, na2, natom, 00627 nb1, nb2, ncoa, ncob, ngpts, nkind, nseta, nsetb, ntasks, nthread, 00628 sgfa, sgfb 00629 INTEGER(kind=int_8), DIMENSION(:), 00630 POINTER :: atom_pair_recv, atom_pair_send 00631 INTEGER(kind=int_8), DIMENSION(:, :), 00632 POINTER :: tasks 00633 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, 00634 lb_min, npgfa, npgfb, nsgfa, 00635 nsgfb 00636 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 00637 LOGICAL :: atom_pair_changed, den_found, den_found_a, 00638 distributed_rs_grids, do_igaim, failure, map_consistent, 00639 my_minimum_image, my_soft 00640 REAL(dp), DIMENSION(:, :, :), POINTER :: my_current, my_gauge, my_rho 00641 REAL(KIND=dp) :: eps_rho_rspace, 00642 kind_radius_a, kind_radius_b, 00643 Lxo2, Lyo2, Lzo2, rab2, 00644 scale, scale2, zetp 00645 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb 00646 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 00647 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dist_ab, jp_block_a, 00648 jp_block_b, jp_block_c, jp_block_d, jpab_a, jpab_b, jpab_c, jpab_d, 00649 jpblock_a, jpblock_b, jpblock_c, jpblock_d, rpgfa, rpgfb, sphi_a, 00650 sphi_b, work, zeta, zetb 00651 REAL(KIND=dp), DIMENSION(:, :, :), 00652 POINTER :: jpabt_a, jpabt_b, jpabt_c, 00653 jpabt_d, workt 00654 TYPE(atomic_kind_type), DIMENSION(:), 00655 POINTER :: atomic_kind_set 00656 TYPE(atomic_kind_type), POINTER :: atomic_kind 00657 TYPE(cell_type), POINTER :: cell 00658 TYPE(cp_dbcsr_type), POINTER :: deltajp_a, deltajp_b, 00659 deltajp_c, deltajp_d, mat_a, 00660 mat_b, mat_c, mat_d 00661 TYPE(cp_para_env_type), POINTER :: para_env 00662 TYPE(cube_info_type), DIMENSION(:), 00663 POINTER :: cube_info 00664 TYPE(dft_control_type), POINTER :: dft_control 00665 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 00666 TYPE(gto_basis_set_p_type), 00667 DIMENSION(:), POINTER :: basis_set_list 00668 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, 00669 orb_basis_set 00670 TYPE(lgrid_type) :: lgrid 00671 TYPE(neighbor_list_iterator_p_type), 00672 DIMENSION(:), POINTER :: nl_iterator 00673 TYPE(neighbor_list_set_p_type), 00674 DIMENSION(:), POINTER :: sab_orb 00675 TYPE(particle_type), DIMENSION(:), 00676 POINTER :: particle_set 00677 TYPE(pw_env_type), POINTER :: pw_env 00678 TYPE(realspace_grid_desc_p_type), 00679 DIMENSION(:), POINTER :: rs_descs 00680 TYPE(realspace_grid_p_type), 00681 DIMENSION(:), POINTER :: rs_current, rs_rho 00682 TYPE(realspaces_grid_p_type), 00683 DIMENSION(:), POINTER :: rs_gauge 00684 TYPE(section_vals_type), POINTER :: input, interp_section 00685 00686 failure=.FALSE. 00687 NULLIFY(atomic_kind, cell, dft_control, orb_basis_set, rs_rho, & 00688 atomic_kind_set, sab_orb, particle_set, rs_current, pw_env, & 00689 rs_descs, para_env, dist_ab, set_radius_a, set_radius_b, la_max, & 00690 la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, rpgfa, rpgfb, & 00691 sphi_a, sphi_b, zeta, zetb, first_sgfa, first_sgfb, dist_ab, & 00692 tasks, workt, mat_a, mat_b, mat_c, mat_d, rs_gauge) 00693 NULLIFY(deltajp_a,deltajp_b,deltajp_c,deltajp_d) 00694 NULLIFY(jp_block_a,jp_block_b,jp_block_c,jp_block_d) 00695 NULLIFY(jpblock_a,jpblock_b,jpblock_c,jpblock_d) 00696 NULLIFY(jpabt_a,jpabt_b,jpabt_c,jpabt_d) 00697 NULLIFY(lgrid%r) 00698 00699 CALL timeset(routineN,handle) 00700 00701 ! 00702 ! Set pointers for the different gauge 00703 do_igaim = current_env%gauge.EQ.current_gauge_atom 00704 00705 IF(do_igaim) THEN 00706 mat_a => mat_jp 00707 mat_b => mat_jp_rii 00708 mat_c => mat_jp_riii 00709 mat_d => mat_d0 00710 ELSE 00711 mat_a => mat_jp 00712 mat_b => mat_jp_rii 00713 mat_c => mat_jp_riii 00714 ENDIF 00715 00716 CALL get_qs_env(qs_env=qs_env,& 00717 atomic_kind_set=atomic_kind_set,& 00718 cell=cell,& 00719 dft_control=dft_control,& 00720 particle_set=particle_set,& 00721 sab_all=sab_orb,& 00722 para_env=para_env,& 00723 input=input,& 00724 pw_env=pw_env,error=error) 00725 00726 CALL get_current_env(current_env=current_env,& 00727 & rs_gauge=rs_gauge,& 00728 & error=error) 00729 00730 ! Component of appearing in the vector product rxp, iiB and iiiB 00731 CALL set_vecp(iB,iiB,iiiB) 00732 ! 00733 ! 00734 scale2 = 0.0_dp 00735 idir2 = 1 00736 IF(idir.NE.iB) THEN 00737 CALL set_vecp_rev(idir,iB,idir2) 00738 scale2 = fac_vecp(idir,iB,idir2) 00739 ENDIF 00740 ! 00741 ! *** assign from pw_env 00742 gridlevel_info=>pw_env%gridlevel_info 00743 cube_info=>pw_env%cube_info 00744 00745 interp_section => section_vals_get_subs_vals(input,"DFT%MGRID%INTERPOLATOR",& 00746 error=error) 00747 ! CALL section_vals_val_get(interp_section,"KIND",i_val=interp_kind,error=error) 00748 00749 ! Check that the neighbor list with all the pairs is associated 00750 CPPrecondition(ASSOCIATED(sab_orb),cp_failure_level,routineP,error,failure) 00751 ! *** set up the pw multi-grids 00752 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 00753 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho, error=error) 00754 00755 distributed_rs_grids=.FALSE. 00756 DO igrid_level=1,gridlevel_info%ngrid_levels 00757 IF(.NOT.ALL(rs_descs(igrid_level)%rs_desc%perd == 1)) THEN 00758 distributed_rs_grids = .TRUE. 00759 ENDIF 00760 ENDDO 00761 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 00762 map_consistent = dft_control%qs_control%map_consistent 00763 nthread = 1 00764 00765 ! *** Allocate work storage *** 00766 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00767 maxco=maxco,& 00768 maxsgf=maxsgf,& 00769 maxsgf_set=maxsgf_set) 00770 my_minimum_image = .TRUE. 00771 ! IF(PRESENT(minimum_image)) THEN 00772 ! my_minimum_image=minimum_image 00773 Lxo2 = SQRT ( SUM ( cell % hmat ( :, 1 ) ** 2 ) )/2.0_dp 00774 Lyo2 = SQRT ( SUM ( cell % hmat ( :, 2 ) ** 2 ) )/2.0_dp 00775 Lzo2 = SQRT ( SUM ( cell % hmat ( :, 3 ) ** 2 ) )/2.0_dp 00776 ! END IF 00777 00778 my_soft=.FALSE. 00779 IF (PRESENT(soft_valid)) my_soft = soft_valid 00780 00781 nkind = SIZE(atomic_kind_set) 00782 00783 CALL reallocate(jpabt_a,1,maxco,1,maxco,0,nthread-1) 00784 CALL reallocate(jpabt_b,1,maxco,1,maxco,0,nthread-1) 00785 CALL reallocate(jpabt_c,1,maxco,1,maxco,0,nthread-1) 00786 CALL reallocate(jpabt_d,1,maxco,1,maxco,0,nthread-1) 00787 CALL reallocate(workt,1,maxco,1,maxsgf_set,0,nthread-1) 00788 CALL reallocate(tasks,1,6,1,max_tasks) 00789 CALL reallocate(dist_ab,1,3,1,max_tasks) 00790 00791 tasks = 0 00792 ntasks = 0 00793 curr_tasks = SIZE(tasks,2) 00794 00795 ! get maximum numbers 00796 natom = SIZE( particle_set ) 00797 maxset=0 00798 maxpgf=0 00799 00800 DO ikind=1,nkind 00801 atomic_kind => atomic_kind_set(ikind) 00802 00803 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00804 orb_basis_set=orb_basis_set) 00805 00806 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 00807 00808 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 00809 npgf=npgfa, nset=nseta ) 00810 00811 maxset=MAX(nseta,maxset) 00812 maxpgf=MAX(MAXVAL(npgfa),maxpgf) 00813 END DO 00814 00815 ! *** Initialize working density matrix *** 00816 00817 ! distributed rs grids require a matrix that will be changed (distribute_tasks) 00818 ! whereas this is not the case for replicated grids 00819 IF (distributed_rs_grids) THEN 00820 ALLOCATE(deltajp_a,deltajp_b,deltajp_c) 00821 CALL cp_dbcsr_init(deltajp_a, error=error) 00822 CALL cp_dbcsr_init(deltajp_b, error=error) 00823 CALL cp_dbcsr_init(deltajp_c, error=error) 00824 IF(do_igaim) THEN 00825 ALLOCATE(deltajp_d) 00826 CALL cp_dbcsr_init(deltajp_d, error=error) 00827 ENDIF 00828 00829 CALL cp_dbcsr_create(deltajp_a, ' deltajp_a ', & 00830 cp_dbcsr_distribution (mat_a), cp_dbcsr_get_matrix_type (mat_a),& 00831 cp_dbcsr_row_block_sizes(mat_a),& 00832 cp_dbcsr_col_block_sizes(mat_a), cp_dbcsr_get_num_blocks(mat_a), cp_dbcsr_get_data_size(mat_a),& 00833 error=error) 00834 CALL cp_dbcsr_create(deltajp_b, ' deltajp_b ', & 00835 cp_dbcsr_distribution (mat_a), cp_dbcsr_get_matrix_type (mat_a),& 00836 cp_dbcsr_row_block_sizes(mat_a),& 00837 cp_dbcsr_col_block_sizes(mat_a), cp_dbcsr_get_num_blocks(mat_a), cp_dbcsr_get_data_size(mat_a),& 00838 error=error) 00839 CALL cp_dbcsr_create(deltajp_c, ' deltajp_c ', & 00840 cp_dbcsr_distribution (mat_a), cp_dbcsr_get_matrix_type (mat_a),& 00841 cp_dbcsr_row_block_sizes(mat_a),& 00842 cp_dbcsr_col_block_sizes(mat_a), cp_dbcsr_get_num_blocks(mat_a), cp_dbcsr_get_data_size(mat_a),& 00843 error=error) 00844 IF(do_igaim) CALL cp_dbcsr_create(deltajp_d, ' deltajp_d ', & 00845 cp_dbcsr_distribution (mat_a), cp_dbcsr_get_matrix_type (mat_a),& 00846 cp_dbcsr_row_block_sizes(mat_a),& 00847 cp_dbcsr_col_block_sizes(mat_a), cp_dbcsr_get_num_blocks(mat_a), cp_dbcsr_get_data_size(mat_a),& 00848 error=error) 00849 ELSE 00850 deltajp_a => mat_a !mat_jp 00851 deltajp_b => mat_b !mat_jp_rii 00852 deltajp_c => mat_c !mat_jp_riii 00853 IF(do_igaim) deltajp_d => mat_d !mat_d0 00854 ENDIF 00855 00856 ALLOCATE (basis_set_list(nkind),STAT=istat) 00857 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00858 DO ikind=1,nkind 00859 atomic_kind => atomic_kind_set(ikind) 00860 CALL get_atomic_kind(atomic_kind=atomic_kind,softb=my_soft,orb_basis_set=basis_set_a) 00861 IF (ASSOCIATED(basis_set_a)) THEN 00862 basis_set_list(ikind)%gto_basis_set => basis_set_a 00863 ELSE 00864 NULLIFY(basis_set_list(ikind)%gto_basis_set) 00865 END IF 00866 END DO 00867 CALL neighbor_list_iterator_create(nl_iterator,sab_orb) 00868 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 00869 CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,iatom=iatom,jatom=jatom,r=rab) 00870 basis_set_a => basis_set_list(ikind)%gto_basis_set 00871 IF (.NOT.ASSOCIATED(basis_set_a)) CYCLE 00872 basis_set_b => basis_set_list(jkind)%gto_basis_set 00873 IF (.NOT.ASSOCIATED(basis_set_b)) CYCLE 00874 ra(:) = pbc(particle_set(iatom)%r,cell) 00875 ! basis ikind 00876 first_sgfa => basis_set_a%first_sgf 00877 la_max => basis_set_a%lmax 00878 la_min => basis_set_a%lmin 00879 npgfa => basis_set_a%npgf 00880 nseta = basis_set_a%nset 00881 nsgfa => basis_set_a%nsgf_set 00882 rpgfa => basis_set_a%pgf_radius 00883 set_radius_a => basis_set_a%set_radius 00884 kind_radius_a= basis_set_a%kind_radius 00885 sphi_a => basis_set_a%sphi 00886 zeta => basis_set_a%zet 00887 ! basis jkind 00888 first_sgfb => basis_set_b%first_sgf 00889 lb_max => basis_set_b%lmax 00890 lb_min => basis_set_b%lmin 00891 npgfb => basis_set_b%npgf 00892 nsetb = basis_set_b%nset 00893 nsgfb => basis_set_b%nsgf_set 00894 rpgfb => basis_set_b%pgf_radius 00895 set_radius_b => basis_set_b%set_radius 00896 kind_radius_b= basis_set_b%kind_radius 00897 sphi_b => basis_set_b%sphi 00898 zetb => basis_set_b%zet 00899 00900 IF(my_minimum_image) THEN 00901 IF(ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) THEN 00902 CYCLE 00903 END IF 00904 END IF 00905 00906 brow = iatom 00907 bcol = jatom 00908 00909 CALL cp_dbcsr_get_block_p(matrix=mat_a,row=brow,col=bcol,& 00910 block=jp_block_a,found=den_found_a) 00911 CALL cp_dbcsr_get_block_p(matrix=mat_b,row=brow,col=bcol,& 00912 block=jp_block_b,found=den_found) 00913 CALL cp_dbcsr_get_block_p(matrix=mat_c,row=brow,col=bcol,& 00914 block=jp_block_c,found=den_found) 00915 IF(do_igaim) CALL cp_dbcsr_get_block_p(matrix=mat_d,row=brow,col=bcol,& 00916 block=jp_block_d,found=den_found) 00917 00918 IF (.NOT.ASSOCIATED(jp_block_a)) CYCLE 00919 00920 IF (distributed_rs_grids) THEN 00921 NULLIFY (jpblock_a,jpblock_b,jpblock_c,jpblock_d) 00922 CALL cp_dbcsr_add_block_node ( deltajp_a, brow, bcol, jpblock_a ,error=error) 00923 jpblock_a = jp_block_a 00924 CALL cp_dbcsr_add_block_node ( deltajp_b, brow, bcol, jpblock_b ,error=error) 00925 jpblock_b = jp_block_b 00926 CALL cp_dbcsr_add_block_node ( deltajp_c, brow, bcol, jpblock_c ,error=error) 00927 jpblock_c = jp_block_c 00928 IF(do_igaim) CALL cp_dbcsr_add_block_node ( deltajp_d, brow, bcol, jpblock_d ,error=error) 00929 IF(do_igaim) jpblock_d = jp_block_d 00930 ELSE 00931 jpblock_a => jp_block_a 00932 jpblock_b => jp_block_b 00933 jpblock_c => jp_block_c 00934 IF(do_igaim) jpblock_d => jp_block_d 00935 ENDIF 00936 00937 IF (.NOT. map_consistent) THEN 00938 IF ( ALL ( 100.0_dp*ABS(jpblock_a) < eps_rho_rspace ) .AND. & 00939 ALL ( 100.0_dp*ABS(jpblock_b) < eps_rho_rspace ) .AND. & 00940 ALL ( 100.0_dp*ABS(jpblock_c) < eps_rho_rspace ) ) THEN 00941 CYCLE 00942 END IF 00943 END IF 00944 00945 CALL task_list_inner_loop(tasks, dist_ab, ntasks, curr_tasks, rs_descs,dft_control,cube_info,gridlevel_info,cell,& 00946 iatom,jatom,rpgfa,rpgfb,zeta,zetb,kind_radius_b,set_radius_a,set_radius_b,ra,rab,& 00947 la_max,la_min,lb_max,lb_min,npgfa,npgfb,maxpgf,maxset,natom,nseta,nsetb,error) 00948 00949 END DO 00950 CALL neighbor_list_iterator_release(nl_iterator) 00951 00952 DEALLOCATE (basis_set_list,STAT=istat) 00953 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00954 00955 IF (distributed_rs_grids) THEN 00956 CALL cp_dbcsr_finalize(deltajp_a, error=error) 00957 CALL cp_dbcsr_finalize(deltajp_b, error=error) 00958 CALL cp_dbcsr_finalize(deltajp_c, error=error) 00959 IF(do_igaim) CALL cp_dbcsr_finalize(deltajp_d, error=error) 00960 ENDIF 00961 00962 ! sorts / redistributes the task list 00963 CALL distribute_tasks ( rs_descs, ntasks, natom,& 00964 maxset, maxpgf, tasks, dist_ab, atom_pair_send, atom_pair_recv,& 00965 symmetric=.FALSE., reorder_rs_grid_ranks=.TRUE., & 00966 skip_load_balance_distributed=.FALSE., error=error) 00967 00968 ALLOCATE(rs_current(gridlevel_info%ngrid_levels)) 00969 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00970 00971 DO igrid_level=1,gridlevel_info%ngrid_levels 00972 ! Here we need to reallocate the distributed rs_grids, which may have been reordered 00973 ! by distribute_tasks 00974 IF (rs_descs(igrid_level)%rs_desc%distributed) THEN 00975 CALL rs_grid_release(rs_rho(igrid_level)%rs_grid, error=error) 00976 NULLIFY (rs_rho(igrid_level)%rs_grid) 00977 CALL rs_grid_create(rs_rho(igrid_level)%rs_grid, rs_descs(igrid_level)%rs_desc, error=error) 00978 ELSE 00979 CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid, error=error) 00980 ENDIF 00981 CALL rs_grid_zero(rs_rho(igrid_level)%rs_grid) 00982 CALL rs_grid_create(rs_current(igrid_level)%rs_grid, rs_descs(igrid_level)%rs_desc,& 00983 & error=error) 00984 CALL rs_grid_zero(rs_current(igrid_level)%rs_grid) 00985 ENDDO 00986 00987 ngpts = 0 00988 IF ( nthread > 1 ) THEN 00989 DO igrid_level = 1,gridlevel_info%ngrid_levels 00990 ngpts = MAX(ngpts,rs_rho(igrid_level)%rs_grid%ngpts_local) 00991 END DO 00992 ngpts = ngpts*nthread 00993 CALL reallocate(lgrid%r,1,ngpts) 00994 END IF 00995 00996 ! 00997 ! we need to build the gauge here 00998 IF(.NOT.current_env%gauge_init) THEN 00999 CALL current_set_gauge(current_env, qs_env, error=error) 01000 current_env%gauge_init = .TRUE. 01001 ENDIF 01002 ! 01003 ! for any case double check the bounds ! 01004 DO igrid_level=1,gridlevel_info%ngrid_levels 01005 my_rho=>rs_rho(igrid_level)%rs_grid%r 01006 my_current=>rs_current(igrid_level)%rs_grid%r 01007 IF(LBOUND( my_rho, 3 ).NE.LBOUND( my_current, 3 ).OR.& 01008 LBOUND( my_rho, 2 ).NE.LBOUND( my_current, 2 ).OR.& 01009 LBOUND( my_rho, 1 ).NE.LBOUND( my_current, 1 ).OR.& 01010 UBOUND( my_rho, 3 ).NE.UBOUND( my_current, 3 ).OR.& 01011 UBOUND( my_rho, 2 ).NE.UBOUND( my_current, 2 ).OR.& 01012 UBOUND( my_rho, 1 ).NE.UBOUND( my_current, 1 ) )THEN 01013 WRITE(*,*) 'LBOUND(my_rho,3),LBOUND(my_current,3)',LBOUND(my_rho,3),LBOUND(my_current,3) 01014 WRITE(*,*) 'LBOUND(my_rho,2),LBOUND(my_current,2)',LBOUND(my_rho,2),LBOUND(my_current,2) 01015 WRITE(*,*) 'LBOUND(my_rho,1),LBOUND(my_current,1)',LBOUND(my_rho,1),LBOUND(my_current,1) 01016 WRITE(*,*) 'UBOUND(my_rho,3),UBOUND(my_current,3)',UBOUND(my_rho,3),UBOUND(my_current,3) 01017 WRITE(*,*) 'UBOUND(my_rho,2),UBOUND(my_current,2)',UBOUND(my_rho,2),UBOUND(my_current,2) 01018 WRITE(*,*) 'UBOUND(my_rho,1),UBOUND(my_current,1)',UBOUND(my_rho,1),UBOUND(my_current,1) 01019 CALL stop_program(routineN,moduleN,__LINE__,"Bug") 01020 ENDIF 01021 01022 IF(do_igaim) THEN 01023 my_gauge=>rs_gauge(1)%rs(igrid_level)%rs_grid%r 01024 IF(LBOUND( my_rho, 3 ).NE.LBOUND( my_gauge, 3 ).OR.& 01025 LBOUND( my_rho, 2 ).NE.LBOUND( my_gauge, 2 ).OR.& 01026 LBOUND( my_rho, 1 ).NE.LBOUND( my_gauge, 1 ).OR.& 01027 UBOUND( my_rho, 3 ).NE.UBOUND( my_gauge, 3 ).OR.& 01028 UBOUND( my_rho, 2 ).NE.UBOUND( my_gauge, 2 ).OR.& 01029 UBOUND( my_rho, 1 ).NE.UBOUND( my_gauge, 1 )) THEN 01030 WRITE(*,*) 'LBOUND(my_rho,3),LBOUND(my_gauge,3)',LBOUND(my_rho,3),LBOUND(my_gauge,3) 01031 WRITE(*,*) 'LBOUND(my_rho,2),LBOUND(my_gauge,2)',LBOUND(my_rho,2),LBOUND(my_gauge,2) 01032 WRITE(*,*) 'LBOUND(my_rho,1),LBOUND(my_gauge,1)',LBOUND(my_rho,1),LBOUND(my_gauge,1) 01033 WRITE(*,*) 'UBOUND(my_rho,3),UbOUND(my_gauge,3)',UBOUND(my_rho,3),UBOUND(my_gauge,3) 01034 WRITE(*,*) 'UBOUND(my_rho,2),UBOUND(my_gauge,2)',UBOUND(my_rho,2),UBOUND(my_gauge,2) 01035 WRITE(*,*) 'UBOUND(my_rho,1),UBOUND(my_gauge,1)',UBOUND(my_rho,1),UBOUND(my_gauge,1) 01036 CALL stop_program(routineN,moduleN,__LINE__,"Bug") 01037 ENDIF 01038 ENDIF 01039 ENDDO 01040 ! 01041 !------------------------------------------------------------- 01042 01043 IF (distributed_rs_grids) THEN 01044 CALL rs_distribute_matrix (rs_descs, deltajp_a, atom_pair_send, atom_pair_recv, & 01045 natom, scatter=.TRUE., error=error) 01046 CALL rs_distribute_matrix (rs_descs, deltajp_b, atom_pair_send, atom_pair_recv, & 01047 natom, scatter=.TRUE., error=error) 01048 CALL rs_distribute_matrix (rs_descs, deltajp_c, atom_pair_send, atom_pair_recv, & 01049 natom, scatter=.TRUE., error=error) 01050 IF(do_igaim) CALL rs_distribute_matrix (rs_descs, deltajp_d, atom_pair_send, atom_pair_recv, & 01051 natom, scatter=.TRUE., error=error) 01052 ENDIF 01053 01054 ithread = 0 01055 jpab_a => jpabt_a(:,:,ithread) 01056 jpab_b => jpabt_b(:,:,ithread) 01057 jpab_c => jpabt_c(:,:,ithread) 01058 IF(do_igaim) jpab_d => jpabt_d(:,:,ithread) 01059 work => workt(:,:,ithread) 01060 01061 iatom_old = -1 ; jatom_old = -1 ; iset_old = -1 ; jset_old = -1 01062 ikind_old = -1 ; jkind_old = -1 01063 01064 loop_tasks: DO itask = 1,ntasks 01065 01066 CALL int2pair(tasks(3,itask),igrid_level,iatom,jatom,iset,jset,ipgf,jpgf,natom,maxset,maxpgf) 01067 01068 ! apparently generalised collocation not implemented correctly yet 01069 CPPostcondition(tasks(4,itask).NE.2,cp_failure_level,routineP,error,failure) 01070 01071 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN 01072 01073 ikind = particle_set(iatom)%atomic_kind%kind_number 01074 jkind = particle_set(jatom)%atomic_kind%kind_number 01075 01076 IF(iatom .NE. iatom_old ) ra(:) = pbc(particle_set(iatom)%r,cell) 01077 01078 brow = iatom 01079 bcol = jatom 01080 01081 IF (ikind .NE. ikind_old ) THEN 01082 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01083 softb = my_soft, & 01084 orb_basis_set=orb_basis_set) 01085 01086 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01087 first_sgf=first_sgfa,& 01088 lmax=la_max,& 01089 lmin=la_min,& 01090 npgf=npgfa,& 01091 nset=nseta,& 01092 nsgf_set=nsgfa,& 01093 pgf_radius=rpgfa,& 01094 set_radius=set_radius_a,& 01095 sphi=sphi_a,& 01096 zet=zeta) 01097 ENDIF 01098 01099 IF (jkind .NE. jkind_old ) THEN 01100 01101 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind,& 01102 softb = my_soft, & 01103 orb_basis_set=orb_basis_set) 01104 01105 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01106 first_sgf=first_sgfb,& 01107 kind_radius=kind_radius_b,& 01108 lmax=lb_max,& 01109 lmin=lb_min,& 01110 npgf=npgfb,& 01111 nset=nsetb,& 01112 nsgf_set=nsgfb,& 01113 pgf_radius=rpgfb,& 01114 set_radius=set_radius_b,& 01115 sphi=sphi_b,& 01116 zet=zetb) 01117 01118 ENDIF 01119 01120 CALL cp_dbcsr_get_block_p(matrix=deltajp_a,row=brow,col=bcol,& 01121 block=jp_block_a,found=den_found) 01122 CALL cp_dbcsr_get_block_p(matrix=deltajp_b,row=brow,col=bcol,& 01123 block=jp_block_b,found=den_found) 01124 CALL cp_dbcsr_get_block_p(matrix=deltajp_c,row=brow,col=bcol,& 01125 block=jp_block_c,found=den_found) 01126 IF(do_igaim) CALL cp_dbcsr_get_block_p(matrix=deltajp_d,row=brow,col=bcol,& 01127 block=jp_block_d,found=den_found) 01128 01129 IF (.NOT.ASSOCIATED(jp_block_a)) & 01130 CALL stop_program(routineN,moduleN,__LINE__,& 01131 "p_block not associated in deltap") 01132 01133 iatom_old = iatom 01134 jatom_old = jatom 01135 ikind_old = ikind 01136 jkind_old = jkind 01137 01138 atom_pair_changed = .TRUE. 01139 01140 ELSE 01141 01142 atom_pair_changed = .FALSE. 01143 01144 ENDIF 01145 01146 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN 01147 01148 ncoa = npgfa(iset)*ncoset(la_max(iset)) 01149 sgfa = first_sgfa(1,iset) 01150 ncob = npgfb(jset)*ncoset(lb_max(jset)) 01151 sgfb = first_sgfb(1,jset) 01152 ! Decontraction step for the selected blocks of the 3 density matrices 01153 01154 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 01155 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01156 jp_block_a(sgfa,sgfb),SIZE(jp_block_a,1),& 01157 0.0_dp,work(1,1),maxco) 01158 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 01159 1.0_dp,work(1,1),maxco,& 01160 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01161 0.0_dp,jpab_a(1,1),maxco) 01162 01163 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 01164 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01165 jp_block_b(sgfa,sgfb),SIZE(jp_block_b,1),& 01166 0.0_dp,work(1,1),maxco) 01167 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 01168 1.0_dp,work(1,1),maxco,& 01169 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01170 0.0_dp,jpab_b(1,1),maxco) 01171 01172 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 01173 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01174 jp_block_c(sgfa,sgfb),SIZE(jp_block_c,1),& 01175 0.0_dp,work(1,1),maxco) 01176 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 01177 1.0_dp,work(1,1),maxco,& 01178 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01179 0.0_dp,jpab_c(1,1),maxco) 01180 01181 IF(do_igaim) THEN 01182 CALL dgemm("N","N",ncoa,nsgfb(jset),nsgfa(iset),& 01183 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01184 jp_block_d(sgfa,sgfb),SIZE(jp_block_d,1),& 01185 0.0_dp,work(1,1),maxco) 01186 CALL dgemm("N","T",ncoa,ncob,nsgfb(jset),& 01187 1.0_dp,work(1,1),maxco,& 01188 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01189 0.0_dp,jpab_d(1,1),maxco) 01190 ENDIF 01191 01192 iset_old = iset 01193 jset_old = jset 01194 01195 ENDIF 01196 01197 rab(:) = dist_ab (:,itask) 01198 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 01199 rb(:) = ra(:) + rab(:) 01200 zetp = zeta(ipgf,iset) + zetb(jpgf,jset) 01201 01202 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 01203 na2 = ipgf*ncoset(la_max(iset)) 01204 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1 01205 nb2 = jpgf*ncoset(lb_max(jset)) 01206 01207 ! Four calls to the general collocate density, to multply the correct function 01208 ! to each density matrix 01209 IF(do_igaim) THEN 01210 ! 01211 ! here the decontracted mat_jp_{ab} is multiplied by 01212 ! f_{ab} = g_{a} (dg_{b}/dr)_{idir} - (dg_{a}/dr)_{idir} g_{b} 01213 scale = 1.0_dp 01214 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01215 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01216 ra,rab,rab2,scale,jpab_a,na1-1,nb1-1,& 01217 rs_current(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01218 eps_rho_rspace,& 01219 ga_gb_function=FUNC_ADBmDAB,& 01220 idir=idir,& 01221 map_consistent=map_consistent,error=error) 01222 ! here the decontracted mat_jb_{ab} is multiplied by 01223 ! f_{ab} = g_{a} * g_{b} ! THIS GOES OUTSIDE THE LOOP ! 01224 IF(scale2.NE.0.0_dp) THEN 01225 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01226 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01227 ra,rab,rab2,scale2,jpab_d,na1-1,nb1-1,& 01228 rs_rho(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01229 eps_rho_rspace,& 01230 ga_gb_function=FUNC_AB,& 01231 map_consistent=map_consistent,error=error) 01232 ENDIF!rm 01233 ! here the decontracted mat_jp_rii{ab} is multiplied by 01234 ! f_{ab} = g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} - 01235 ! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b} 01236 scale = 1.0_dp 01237 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01238 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01239 ra,rab,rab2,scale,jpab_b,na1-1,nb1-1,& 01240 rs_current(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01241 eps_rho_rspace,& 01242 ga_gb_function=FUNC_ADBmDAB,& 01243 idir=idir,ir=iiiB,& 01244 rsgauge=rs_gauge(iiiB)%rs(igrid_level)%rs_grid,& 01245 rsbuf=current_env%rs_buf(igrid_level)%rs_grid,& 01246 map_consistent=map_consistent,error=error) 01247 ! here the decontracted mat_jp_riii{ab} is multiplied by 01248 ! f_{ab} = -g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} + 01249 ! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b} 01250 scale = -1.0_dp 01251 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01252 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01253 ra,rab,rab2,scale,jpab_c,na1-1,nb1-1,& 01254 rs_current(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01255 eps_rho_rspace,& 01256 ga_gb_function=FUNC_ADBmDAB,& 01257 idir=idir,ir=iiB,& 01258 rsgauge=rs_gauge(iiB)%rs(igrid_level)%rs_grid,& 01259 rsbuf=current_env%rs_buf(igrid_level)%rs_grid,& 01260 map_consistent=map_consistent,error=error) 01261 ELSE 01262 ! here the decontracted mat_jp_{ab} is multiplied by 01263 ! f_{ab} = g_{a} (dg_{b}/dr)_{idir} - (dg_{a}/dr)_{idir} g_{b} 01264 scale = 1.0_dp 01265 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01266 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01267 ra,rab,rab2,scale,jpab_a,na1-1,nb1-1,& 01268 rs_current(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01269 eps_rho_rspace,& 01270 ga_gb_function=FUNC_ADBmDAB,& 01271 idir=idir,& 01272 map_consistent=map_consistent,error=error) 01273 ! here the decontracted mat_jp_rii{ab} is multiplied by 01274 ! f_{ab} = g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} - 01275 ! (dg_{a}/dr)_{idir} (r - R_{b})_{iiB} g_{b} 01276 scale = 1.0_dp 01277 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01278 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01279 ra,rab,rab2,scale,jpab_b,na1-1,nb1-1,& 01280 rs_current(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01281 eps_rho_rspace,& 01282 ga_gb_function=FUNC_ARDBmDARB,& 01283 idir=idir,ir=iiiB,& 01284 map_consistent=map_consistent,error=error) 01285 ! here the decontracted mat_jp_riii{ab} is multiplied by 01286 ! f_{ab} = -g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} + 01287 ! (dg_{a}/dr)_{idir} (r - R_{b})_{iiB} g_{b} 01288 scale = -1.0_dp 01289 CALL collocate_pgf_product_rspace(la_max(iset),zeta(ipgf,iset),& 01290 la_min(iset),lb_max(jset),zetb(jpgf,jset),lb_min(jset),& 01291 ra,rab,rab2,scale,jpab_c,na1-1,nb1-1,& 01292 rs_current(igrid_level)%rs_grid,cell,cube_info(igrid_level),& 01293 eps_rho_rspace,& 01294 ga_gb_function=FUNC_ARDBmDARB,& 01295 idir=idir,ir=iiB,& 01296 map_consistent=map_consistent,error=error) 01297 ENDIF 01298 01299 END DO loop_tasks 01300 ! 01301 ! Scale the density with the gauge rho * ( r - d(r) ) if needed 01302 IF(do_igaim) THEN 01303 DO igrid_level=1,gridlevel_info%ngrid_levels 01304 CALL rs_grid_mult_and_add(rs_current(igrid_level)%rs_grid,rs_rho(igrid_level)%rs_grid,& 01305 & rs_gauge(idir2)%rs(igrid_level)%rs_grid,1.0_dp) 01306 ENDDO 01307 ENDIF 01308 ! *** Release work storage *** 01309 01310 IF (distributed_rs_grids) THEN 01311 CALL cp_dbcsr_deallocate_matrix ( deltajp_a ,error=error) 01312 CALL cp_dbcsr_deallocate_matrix ( deltajp_b ,error=error) 01313 CALL cp_dbcsr_deallocate_matrix ( deltajp_c ,error=error) 01314 IF(do_igaim) CALL cp_dbcsr_deallocate_matrix ( deltajp_d ,error=error) 01315 END IF 01316 01317 IF ( nthread > 1 ) THEN 01318 DEALLOCATE (lgrid%r,STAT=istat) 01319 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01320 END IF 01321 01322 DEALLOCATE (jpabt_a,jpabt_b,jpabt_c,jpabt_d,workt,tasks,dist_ab,STAT=istat) 01323 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01324 01325 IF (distributed_rs_grids) THEN 01326 DEALLOCATE(atom_pair_send,atom_pair_recv,STAT=istat) 01327 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01328 ENDIF 01329 01330 CALL density_rs2pw(pw_env,rs_current,current_rs,current_gs,& 01331 interp_section=interp_section,error=error) 01332 01333 IF (ASSOCIATED(rs_rho)) THEN 01334 DO i=1, SIZE(rs_rho) 01335 CALL rs_grid_release(rs_rho(i)%rs_grid, error=error) 01336 END DO 01337 END IF 01338 01339 ! Free the array of grids (grids themselves are released in density_rs2pw) 01340 DEALLOCATE (rs_current,STAT=istat) 01341 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01342 01343 CALL timestop(handle) 01344 01345 END SUBROUTINE calculate_jrho_resp 01346 01347 01348 SUBROUTINE current_set_gauge(current_env,qs_env,error) 01349 ! 01350 TYPE(current_env_type) :: current_env 01351 TYPE(qs_environment_type), POINTER :: qs_env 01352 TYPE(cp_error_type), INTENT(INOUT) :: error 01353 01354 CHARACTER(LEN=*), PARAMETER :: routineN = 'current_set_gauge', 01355 routineP = moduleN//':'//routineN 01356 01357 REAL(dp) :: dbox(3) 01358 REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: box_data 01359 INTEGER :: handle, igrid_level, nbox(3), gauge, istat 01360 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: box_ptr 01361 TYPE(realspace_grid_desc_p_type), DIMENSION(:), 01362 POINTER :: rs_descs 01363 TYPE(pw_env_type), POINTER :: pw_env 01364 TYPE(realspaces_grid_p_type), DIMENSION(:), POINTER :: rs_gauge 01365 01366 TYPE(box_type), DIMENSION(:,:,:), POINTER :: box 01367 LOGICAL :: use_old_gauge_atom, failure 01368 01369 NULLIFY(rs_gauge,box) 01370 01371 CALL timeset(routineN,handle) 01372 01373 CALL get_current_env(current_env=current_env,& 01374 & use_old_gauge_atom=use_old_gauge_atom,& 01375 & rs_gauge=rs_gauge,& 01376 & gauge=gauge,& 01377 & error=error) 01378 01379 IF(gauge.EQ.current_gauge_atom) THEN 01380 CALL get_qs_env(qs_env=qs_env,& 01381 pw_env=pw_env,& 01382 error=error) 01383 CALL pw_env_get(pw_env=pw_env,& 01384 rs_descs=rs_descs,& 01385 error=error) 01386 ! 01387 ! box the atoms 01388 IF(use_old_gauge_atom) THEN 01389 CALL box_atoms(qs_env,error) 01390 ELSE 01391 CALL box_atoms_new(current_env,qs_env,box,error) 01392 ENDIF 01393 ! 01394 ! allocate and build the gauge 01395 ALLOCATE (rs_gauge(1)%rs(pw_env%gridlevel_info%ngrid_levels),STAT=istat) 01396 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01397 ALLOCATE (rs_gauge(2)%rs(pw_env%gridlevel_info%ngrid_levels),STAT=istat) 01398 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01399 ALLOCATE (rs_gauge(3)%rs(pw_env%gridlevel_info%ngrid_levels),STAT=istat) 01400 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01401 DO igrid_level = pw_env%gridlevel_info%ngrid_levels,1,-1 01402 01403 CALL rs_grid_create(rs_gauge(1)%rs(igrid_level)%rs_grid, rs_descs(igrid_level)%rs_desc, error=error) 01404 CALL rs_grid_create(rs_gauge(2)%rs(igrid_level)%rs_grid, rs_descs(igrid_level)%rs_desc, error=error) 01405 CALL rs_grid_create(rs_gauge(3)%rs(igrid_level)%rs_grid, rs_descs(igrid_level)%rs_desc, error=error) 01406 01407 IF(use_old_gauge_atom) THEN 01408 CALL collocate_gauge(current_env,qs_env,& 01409 & rs_gauge(1)%rs(igrid_level)%rs_grid,& 01410 & rs_gauge(2)%rs(igrid_level)%rs_grid,& 01411 & rs_gauge(3)%rs(igrid_level)%rs_grid,& 01412 & error) 01413 ELSE 01414 CALL collocate_gauge_new(current_env,qs_env,& 01415 & rs_gauge(1)%rs(igrid_level)%rs_grid,& 01416 & rs_gauge(2)%rs(igrid_level)%rs_grid,& 01417 & rs_gauge(3)%rs(igrid_level)%rs_grid,& 01418 & box,error) 01419 ENDIF 01420 ENDDO 01421 ! 01422 ! allocate the buf 01423 ALLOCATE (current_env%rs_buf(pw_env%gridlevel_info%ngrid_levels),STAT=istat) 01424 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01425 DO igrid_level = 1,pw_env%gridlevel_info%ngrid_levels 01426 CALL rs_grid_create(current_env%rs_buf(igrid_level)%rs_grid, rs_descs(igrid_level)%rs_desc, error=error) 01427 END DO 01428 ! 01429 DEALLOCATE(box_ptr,box_data) 01430 CALL deallocate_box(box) 01431 ENDIF 01432 01433 CALL timestop(handle) 01434 01435 CONTAINS 01436 01437 SUBROUTINE box_atoms(qs_env,error) 01438 TYPE(qs_environment_type), POINTER :: qs_env 01439 TYPE(cp_error_type), INTENT(INOUT) :: error 01440 01441 REAL(kind=dp), PARAMETER :: box_size_guess = 5.0_dp 01442 01443 INTEGER :: i, iatom, ibox, ii, jbox, 01444 kbox, natms 01445 REAL(dp) :: offset(3) 01446 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom 01447 TYPE(atomic_kind_type), DIMENSION(:), 01448 POINTER :: atomic_kind_set 01449 TYPE(cell_type), POINTER :: cell 01450 TYPE(particle_type), DIMENSION(:), 01451 POINTER :: particle_set 01452 01453 CALL get_qs_env(qs_env=qs_env,& 01454 atomic_kind_set=atomic_kind_set,& 01455 cell=cell,& 01456 particle_set=particle_set,& 01457 error=error) 01458 01459 natms = SIZE(particle_set,1) 01460 ALLOCATE(ratom(3,natms)) 01461 ! 01462 ! box the atoms 01463 nbox(1) = CEILING(cell%hmat(1,1) / box_size_guess) 01464 nbox(2) = CEILING(cell%hmat(2,2) / box_size_guess) 01465 nbox(3) = CEILING(cell%hmat(3,3) / box_size_guess) 01466 !write(*,*) 'nbox',nbox 01467 dbox(1) = cell%hmat(1,1) / REAL(nbox(1),dp) 01468 dbox(2) = cell%hmat(2,2) / REAL(nbox(2),dp) 01469 dbox(3) = cell%hmat(3,3) / REAL(nbox(3),dp) 01470 !write(*,*) 'dbox',dbox 01471 ALLOCATE(box_ptr(0:nbox(1),0:nbox(2)-1,0:nbox(3)-1),box_data(3,natms)) 01472 box_data(:,:) = HUGE(0.0_dp) 01473 box_ptr(:,:,:) = HUGE(0) 01474 ! 01475 offset(1) = cell%hmat(1,1)*0.5_dp 01476 offset(2) = cell%hmat(2,2)*0.5_dp 01477 offset(3) = cell%hmat(3,3)*0.5_dp 01478 DO iatom = 1,natms 01479 ratom(:,iatom) = pbc(particle_set(iatom)%r(:),cell)+offset(:) 01480 ENDDO 01481 ! 01482 i = 1 01483 DO kbox = 0,nbox(3)-1 01484 DO jbox = 0,nbox(2)-1 01485 box_ptr(0,jbox,kbox) = i 01486 DO ibox = 0,nbox(1)-1 01487 ii = 0 01488 DO iatom = 1,natms 01489 IF(INT(ratom(1,iatom)/dbox(1)).EQ.ibox.AND.& 01490 INT(ratom(2,iatom)/dbox(2)).EQ.jbox.AND.& 01491 INT(ratom(3,iatom)/dbox(3)).EQ.kbox) THEN 01492 box_data(:,i) = ratom(:,iatom)-offset(:) 01493 i = i + 1 01494 ii = ii + 1 01495 ENDIF 01496 ENDDO 01497 box_ptr(ibox+1,jbox,kbox) = box_ptr(ibox,jbox,kbox) + ii 01498 ENDDO 01499 ENDDO 01500 ENDDO 01501 ! 01502 IF(.FALSE.) THEN 01503 DO kbox = 0,nbox(3)-1 01504 DO jbox = 0,nbox(2)-1 01505 DO ibox = 0,nbox(1)-1 01506 WRITE(*,*) 'box=',ibox,jbox,kbox 01507 WRITE(*,*) 'nbr atom=',box_ptr(ibox+1,jbox,kbox)-box_ptr(ibox,jbox,kbox) 01508 DO iatom = box_ptr(ibox,jbox,kbox),box_ptr(ibox+1,jbox,kbox)-1 01509 WRITE(*,*) 'iatom=',iatom 01510 WRITE(*,'(A,3E14.6)') 'coor=',box_data(:,iatom) 01511 ENDDO 01512 ENDDO 01513 ENDDO 01514 ENDDO 01515 ENDIF 01516 DEALLOCATE(ratom) 01517 END SUBROUTINE box_atoms 01518 01519 01520 SUBROUTINE collocate_gauge(current_env,qs_env,rs_grid_x,rs_grid_y,rs_grid_z,error) 01521 ! 01522 TYPE(current_env_type) :: current_env 01523 TYPE(qs_environment_type), POINTER :: qs_env 01524 TYPE(realspace_grid_type), POINTER :: rs_grid_x, rs_grid_y, 01525 rs_grid_z 01526 TYPE(cp_error_type), INTENT(INOUT) :: error 01527 01528 INTEGER :: i, iatom, ibeg, ibox, iend, imax, imin, j, jatom, jbox, jmax, 01529 jmin, k, kbox, kmax, kmin, lb(3), lb_local(3), natms, natms_local, ng(3) 01530 REAL(KIND=dp) :: ab, buf_tmp, dist, dr(3), gauge_atom_radius, offset(3), 01531 pa, pb, point(3), pra(3), r(3), res(3), summe, tmp, x, y, z 01532 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, nrm_atms_pnt 01533 REAL(kind=dp), ALLOCATABLE, 01534 DIMENSION(:, :) :: atms_pnt, ratom 01535 REAL(kind=dp), DIMENSION(:, :, :), 01536 POINTER :: grid_x, grid_y, grid_z 01537 TYPE(atomic_kind_type), DIMENSION(:), 01538 POINTER :: atomic_kind_set 01539 TYPE(cell_type), POINTER :: cell 01540 TYPE(particle_type), DIMENSION(:), 01541 POINTER :: particle_set 01542 01543 ! 01544 01545 CALL get_current_env(current_env=current_env,& 01546 gauge_atom_radius=gauge_atom_radius,& 01547 error=error) 01548 ! 01549 CALL get_qs_env(qs_env=qs_env,& 01550 atomic_kind_set=atomic_kind_set,& 01551 cell=cell,& 01552 particle_set=particle_set,& 01553 error=error) 01554 ! 01555 natms = SIZE(particle_set,1) 01556 dr(1) = rs_grid_x%desc%dh(1,1) 01557 dr(2) = rs_grid_x%desc%dh(2,2) 01558 dr(3) = rs_grid_x%desc%dh(3,3) 01559 lb(:) = rs_grid_x%desc%lb(:) 01560 lb_local(:) = rs_grid_x%lb_local(:) 01561 grid_x => rs_grid_x%r(:,:,:) 01562 grid_y => rs_grid_y%r(:,:,:) 01563 grid_z => rs_grid_z%r(:,:,:) 01564 ng(:) = UBOUND(grid_x) 01565 offset(1) = cell%hmat(1,1)*0.5_dp 01566 offset(2) = cell%hmat(2,2)*0.5_dp 01567 offset(3) = cell%hmat(3,3)*0.5_dp 01568 ALLOCATE(buf(natms),ratom(3,natms),atms_pnt(3,natms),nrm_atms_pnt(natms)) 01569 ! 01570 ! go over the grid 01571 DO k = 1,ng(3) 01572 DO j = 1,ng(2) 01573 DO i = 1,ng(1) 01574 ! 01575 point(3) = REAL(k-1+lb_local(3)-lb(3),dp)*dr(3) 01576 point(2) = REAL(j-1+lb_local(2)-lb(2),dp)*dr(2) 01577 point(1) = REAL(i-1+lb_local(1)-lb(1),dp)*dr(1) 01578 point = pbc(point,cell) 01579 ! 01580 ! run over the overlaping boxes 01581 natms_local = 0 01582 kmin = INT((point(3)+offset(3)-gauge_atom_radius)/dbox(3)) 01583 kmax = INT((point(3)+offset(3)+gauge_atom_radius)/dbox(3)) 01584 IF(kmax-kmin+1.GT.nbox(3)) THEN 01585 kmin = 0 01586 kmax = nbox(3)-1 01587 ENDIF 01588 DO kbox = kmin,kmax 01589 jmin = INT((point(2)+offset(2)-gauge_atom_radius)/dbox(2)) 01590 jmax = INT((point(2)+offset(2)+gauge_atom_radius)/dbox(2)) 01591 IF(jmax-jmin+1.GT.nbox(2)) THEN 01592 jmin = 0 01593 jmax = nbox(2)-1 01594 ENDIF 01595 DO jbox = jmin,jmax 01596 imin = INT((point(1)+offset(1)-gauge_atom_radius)/dbox(1)) 01597 imax = INT((point(1)+offset(1)+gauge_atom_radius)/dbox(1)) 01598 IF(imax-imin+1.GT.nbox(1)) THEN 01599 imin = 0 01600 imax = nbox(1)-1 01601 ENDIF 01602 DO ibox = imin,imax 01603 ibeg = box_ptr(MODULO(ibox,nbox(1)) ,MODULO(jbox,nbox(2)),MODULO(kbox,nbox(3))) 01604 iend = box_ptr(MODULO(ibox,nbox(1))+1,MODULO(jbox,nbox(2)),MODULO(kbox,nbox(3)))-1 01605 DO iatom = ibeg,iend 01606 r(:) = pbc(box_data(:,iatom)-point(:),cell) + point(:) 01607 dist = (r(1)-point(1))**2 + (r(2)-point(2))**2 + (r(3)-point(3))**2 01608 IF(dist.LT.gauge_atom_radius**2) THEN 01609 natms_local = natms_local + 1 01610 ratom(:,natms_local) = r(:) 01611 ! 01612 ! compute the distance atoms-point 01613 x = point(1) - r(1) 01614 y = point(2) - r(2) 01615 z = point(3) - r(3) 01616 atms_pnt(1,natms_local) = x 01617 atms_pnt(2,natms_local) = y 01618 atms_pnt(3,natms_local) = z 01619 nrm_atms_pnt(natms_local) = SQRT( x*x + y*y + z*z ) 01620 ENDIF 01621 ENDDO 01622 ENDDO 01623 ENDDO 01624 ENDDO 01625 ! 01626 IF(natms_local.GT.0) THEN 01627 ! 01628 ! 01629 DO iatom = 1,natms_local 01630 buf_tmp = 1.0_dp 01631 pra(1) = atms_pnt(1,iatom) 01632 pra(2) = atms_pnt(2,iatom) 01633 pra(3) = atms_pnt(3,iatom) 01634 pa = nrm_atms_pnt(iatom) 01635 DO jatom = 1,natms_local 01636 IF(iatom.EQ.jatom) CYCLE 01637 pb = nrm_atms_pnt(jatom) 01638 x = pra(1) - atms_pnt(1,jatom) 01639 y = pra(2) - atms_pnt(2,jatom) 01640 z = pra(3) - atms_pnt(3,jatom) 01641 ab = SQRT( x*x + y*y + z*z ) 01642 ! 01643 tmp = ( pa - pb ) / ab 01644 tmp = 0.5_dp * ( 3.0_dp - tmp * tmp ) * tmp 01645 tmp = 0.5_dp * ( 3.0_dp - tmp * tmp ) * tmp 01646 tmp = 0.5_dp * ( 3.0_dp - tmp * tmp ) * tmp 01647 buf_tmp = buf_tmp * 0.5_dp * ( 1.0_dp - tmp ) 01648 ENDDO 01649 buf(iatom) = buf_tmp 01650 ENDDO 01651 res(1) = 0.0_dp 01652 res(2) = 0.0_dp 01653 res(3) = 0.0_dp 01654 summe = 0.0_dp 01655 DO iatom = 1,natms_local 01656 res(1) = res(1) + ratom(1,iatom) * buf(iatom) 01657 res(2) = res(2) + ratom(2,iatom) * buf(iatom) 01658 res(3) = res(3) + ratom(3,iatom) * buf(iatom) 01659 summe = summe + buf(iatom) 01660 ENDDO 01661 res(1) = res(1) / summe 01662 res(2) = res(2) / summe 01663 res(3) = res(3) / summe 01664 grid_x(i,j,k) = point(1) - res(1) 01665 grid_y(i,j,k) = point(2) - res(2) 01666 grid_z(i,j,k) = point(3) - res(3) 01667 ELSE 01668 grid_x(i,j,k) = 0.0_dp 01669 grid_y(i,j,k) = 0.0_dp 01670 grid_z(i,j,k) = 0.0_dp 01671 ENDIF 01672 ENDDO 01673 ENDDO 01674 ENDDO 01675 01676 DEALLOCATE(buf,ratom,atms_pnt,nrm_atms_pnt) 01677 01678 END SUBROUTINE collocate_gauge 01679 01680 SUBROUTINE box_atoms_new(current_env,qs_env,box,error) 01681 TYPE(current_env_type) :: current_env 01682 TYPE(qs_environment_type), POINTER :: qs_env 01683 TYPE(box_type), DIMENSION(:, :, :), 01684 POINTER :: box 01685 TYPE(cp_error_type), INTENT(INOUT) :: error 01686 01687 CHARACTER(LEN=*), PARAMETER :: routineN = 'box_atoms_new', 01688 routineP = moduleN//':'//routineN 01689 01690 INTEGER :: handle, i, iatom, ibeg, ibox, iend, ifind, ii, imax, imin, j, 01691 jatom, jbox, jmax, jmin, k, kbox, kmax, kmin, natms, natms_local 01692 REAL(dp) :: gauge_atom_radius, offset(3), 01693 scale 01694 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom 01695 REAL(dp), DIMENSION(:, :), POINTER :: r_ptr 01696 REAL(kind=dp) :: box_center(3), 01697 box_center_wrap(3), 01698 box_size_guess, r(3) 01699 TYPE(atomic_kind_type), DIMENSION(:), 01700 POINTER :: atomic_kind_set 01701 TYPE(cell_type), POINTER :: cell 01702 TYPE(particle_type), DIMENSION(:), 01703 POINTER :: particle_set 01704 01705 CALL timeset(routineN,handle) 01706 01707 CALL get_qs_env(qs_env=qs_env,& 01708 atomic_kind_set=atomic_kind_set,& 01709 cell=cell,& 01710 particle_set=particle_set,& 01711 error=error) 01712 01713 CALL get_current_env(current_env=current_env,& 01714 gauge_atom_radius=gauge_atom_radius,& 01715 error=error) 01716 01717 scale = 2.0_dp 01718 01719 box_size_guess = gauge_atom_radius/scale 01720 01721 natms = SIZE(particle_set,1) 01722 ALLOCATE(ratom(3,natms)) 01723 01724 ! 01725 ! box the atoms 01726 nbox(1) = CEILING(cell%hmat(1,1) / box_size_guess) 01727 nbox(2) = CEILING(cell%hmat(2,2) / box_size_guess) 01728 nbox(3) = CEILING(cell%hmat(3,3) / box_size_guess) 01729 dbox(1) = cell%hmat(1,1) / REAL(nbox(1),dp) 01730 dbox(2) = cell%hmat(2,2) / REAL(nbox(2),dp) 01731 dbox(3) = cell%hmat(3,3) / REAL(nbox(3),dp) 01732 ALLOCATE(box_ptr(0:nbox(1),0:nbox(2)-1,0:nbox(3)-1),box_data(3,natms)) 01733 box_data(:,:) = HUGE(0.0_dp) 01734 box_ptr(:,:,:) = HUGE(0) 01735 ! 01736 offset(1) = cell%hmat(1,1)*0.5_dp 01737 offset(2) = cell%hmat(2,2)*0.5_dp 01738 offset(3) = cell%hmat(3,3)*0.5_dp 01739 DO iatom = 1,natms 01740 ratom(:,iatom) = pbc(particle_set(iatom)%r(:),cell) 01741 ENDDO 01742 ! 01743 i = 1 01744 DO kbox = 0,nbox(3)-1 01745 DO jbox = 0,nbox(2)-1 01746 box_ptr(0,jbox,kbox) = i 01747 DO ibox = 0,nbox(1)-1 01748 ii = 0 01749 DO iatom = 1,natms 01750 IF(MODULO(FLOOR(ratom(1,iatom)/dbox(1)),nbox(1)).EQ.ibox.AND.& 01751 MODULO(FLOOR(ratom(2,iatom)/dbox(2)),nbox(2)).EQ.jbox.AND.& 01752 MODULO(FLOOR(ratom(3,iatom)/dbox(3)),nbox(3)).EQ.kbox) THEN 01753 box_data(:,i) = ratom(:,iatom) 01754 i = i + 1 01755 ii = ii + 1 01756 ENDIF 01757 ENDDO 01758 box_ptr(ibox+1,jbox,kbox) = box_ptr(ibox,jbox,kbox) + ii 01759 ENDDO 01760 ENDDO 01761 ENDDO 01762 ! 01763 IF(.FALSE.) THEN 01764 DO kbox = 0,nbox(3)-1 01765 DO jbox = 0,nbox(2)-1 01766 DO ibox = 0,nbox(1)-1 01767 IF(box_ptr(ibox+1,jbox,kbox)-box_ptr(ibox,jbox,kbox).GT.0) THEN 01768 WRITE(*,*) 'box=',ibox,jbox,kbox 01769 WRITE(*,*) 'nbr atom=',box_ptr(ibox+1,jbox,kbox)-box_ptr(ibox,jbox,kbox) 01770 DO iatom = box_ptr(ibox,jbox,kbox),box_ptr(ibox+1,jbox,kbox)-1 01771 WRITE(*,'(A,I3,3E14.6)') 'coor=',iatom,box_data(:,iatom) 01772 ENDDO 01773 ENDIF 01774 ENDDO 01775 ENDDO 01776 ENDDO 01777 ENDIF 01778 ! 01779 NULLIFY(box) 01780 ALLOCATE(box(0:nbox(1)-1,0:nbox(2)-1,0:nbox(3)-1)) 01781 ! 01782 ! build the list 01783 DO k = 0,nbox(3)-1 01784 DO j = 0,nbox(2)-1 01785 DO i = 0,nbox(1)-1 01786 ! 01787 box_center(1) = ( REAL(i,dp) + 0.5_dp ) * dbox(1) 01788 box_center(2) = ( REAL(j,dp) + 0.5_dp ) * dbox(2) 01789 box_center(3) = ( REAL(k,dp) + 0.5_dp ) * dbox(3) 01790 box_center_wrap = pbc(box_center,cell) 01791 ! 01792 ! find the atoms that are in the overlaping boxes 01793 natms_local = 0 01794 kmin = FLOOR( ( box_center(3) - gauge_atom_radius ) / dbox(3) ) 01795 kmax = FLOOR( ( box_center(3) + gauge_atom_radius ) / dbox(3) ) 01796 IF(kmax-kmin+1.GT.nbox(3)) THEN 01797 kmin = 0 01798 kmax = nbox(3)-1 01799 ENDIF 01800 DO kbox = kmin,kmax 01801 jmin = FLOOR( ( box_center(2) - gauge_atom_radius ) / dbox(2) ) 01802 jmax = FLOOR( ( box_center(2) + gauge_atom_radius ) / dbox(2) ) 01803 IF(jmax-jmin+1.GT.nbox(2)) THEN 01804 jmin = 0 01805 jmax = nbox(2)-1 01806 ENDIF 01807 DO jbox = jmin,jmax 01808 imin = FLOOR( ( box_center(1) - gauge_atom_radius ) / dbox(1) ) 01809 imax = FLOOR( ( box_center(1) + gauge_atom_radius ) / dbox(1) ) 01810 IF(imax-imin+1.GT.nbox(1)) THEN 01811 imin = 0 01812 imax = nbox(1)-1 01813 ENDIF 01814 DO ibox = imin,imax 01815 ibeg = box_ptr(MODULO(ibox,nbox(1)) ,MODULO(jbox,nbox(2)),MODULO(kbox,nbox(3))) 01816 iend = box_ptr(MODULO(ibox,nbox(1))+1,MODULO(jbox,nbox(2)),MODULO(kbox,nbox(3)))-1 01817 DO iatom = ibeg,iend 01818 r = pbc(box_center_wrap(:)-box_data(:,iatom),cell) 01819 IF(ABS(r(1)).LE.(scale+0.5_dp)*dbox(1).AND.& 01820 ABS(r(2)).LE.(scale+0.5_dp)*dbox(2).AND.& 01821 ABS(r(3)).LE.(scale+0.5_dp)*dbox(3)) THEN 01822 natms_local = natms_local + 1 01823 ratom(:,natms_local) = box_data(:,iatom) 01824 ENDIF 01825 ENDDO 01826 ENDDO ! box 01827 ENDDO 01828 ENDDO 01829 ! 01830 ! set the list 01831 box(i,j,k)%n = natms_local 01832 NULLIFY(box(i,j,k)%r) 01833 IF(natms_local.GT.0) THEN 01834 ALLOCATE(box(i,j,k)%r(3,natms_local)) 01835 r_ptr => box(i,j,k)%r 01836 CALL dcopy(3*natms_local,ratom(1,1),1,r_ptr(1,1),1) 01837 ENDIF 01838 ENDDO ! list 01839 ENDDO 01840 ENDDO 01841 01842 IF(.FALSE.) THEN 01843 DO k = 0,nbox(3)-1 01844 DO j = 0,nbox(2)-1 01845 DO i = 0,nbox(1)-1 01846 IF(box(i,j,k)%n.GT.0) THEN 01847 WRITE(*,*) 01848 WRITE(*,*) 'box=',i,j,k 01849 box_center(1) = ( REAL(i,dp)+0.5_dp ) * dbox(1) 01850 box_center(2) = ( REAL(j,dp)+0.5_dp ) * dbox(2) 01851 box_center(3) = ( REAL(k,dp)+0.5_dp ) * dbox(3) 01852 box_center = pbc(box_center,cell) 01853 WRITE(*,'(A,3E14.6)') 'box_center=',box_center 01854 WRITE(*,*) 'nbr atom=',box(i,j,k)%n 01855 r_ptr => box(i,j,k)%r 01856 DO iatom = 1,box(i,j,k)%n 01857 WRITE(*,'(A,I3,3E14.6)') 'coor=',iatom,r_ptr(:,iatom) 01858 r(:) = pbc(box_center(:)-r_ptr(:,iatom),cell) 01859 IF(ABS(r(1)).GT.(scale+0.5_dp)*dbox(1).OR.& 01860 ABS(r(2)).GT.(scale+0.5_dp)*dbox(2).OR.& 01861 ABS(r(3)).GT.(scale+0.5_dp)*dbox(3)) THEN 01862 WRITE(*,*) 'error too many atoms' 01863 WRITE(*,*) 'dist=',ABS(r(:)) 01864 WRITE(*,*) 'large_dist=',(scale+0.5_dp)*dbox 01865 STOP 01866 ENDIF 01867 ENDDO 01868 ENDIF 01869 ENDDO ! list 01870 ENDDO 01871 ENDDO 01872 ENDIF 01873 01874 IF(.TRUE.) THEN 01875 DO k = 0,nbox(3)-1 01876 DO j = 0,nbox(2)-1 01877 DO i = 0,nbox(1)-1 01878 box_center(1) = ( REAL(i,dp)+0.5_dp ) * dbox(1) 01879 box_center(2) = ( REAL(j,dp)+0.5_dp ) * dbox(2) 01880 box_center(3) = ( REAL(k,dp)+0.5_dp ) * dbox(3) 01881 box_center = pbc(box_center,cell) 01882 r_ptr => box(i,j,k)%r 01883 DO iatom = 1,natms 01884 r(:) = pbc(box_center(:)-ratom(:,iatom),cell) 01885 ifind = 0 01886 DO jatom = 1,box(i,j,k)%n 01887 IF(SUM(ABS(ratom(:,iatom)-r_ptr(:,jatom))).LT.1E-10_dp) ifind = 1 01888 ENDDO 01889 01890 IF(ifind.EQ.0) THEN 01891 IF(SQRT(DOT_PRODUCT(r,r)).LT.gauge_atom_radius) THEN 01892 WRITE(*,*) 'error atom too close' 01893 WRITE(*,*) 'iatom',iatom 01894 WRITE(*,*) 'box_center',box_center 01895 WRITE(*,*) 'ratom',ratom(:,iatom) 01896 WRITE(*,*) 'gauge_atom_radius',gauge_atom_radius 01897 STOP 01898 ENDIF 01899 ENDIF 01900 ENDDO 01901 ENDDO ! list 01902 ENDDO 01903 ENDDO 01904 ENDIF 01905 01906 DEALLOCATE(ratom) 01907 01908 CALL timestop(handle) 01909 01910 END SUBROUTINE box_atoms_new 01911 01912 SUBROUTINE collocate_gauge_new(current_env,qs_env,rs_grid_x,rs_grid_y,rs_grid_z,box,error) 01913 ! 01914 TYPE(current_env_type) :: current_env 01915 TYPE(qs_environment_type), POINTER :: qs_env 01916 TYPE(realspace_grid_type), POINTER :: rs_grid_x, rs_grid_y, 01917 rs_grid_z 01918 TYPE(box_type), DIMENSION(:, :, :), 01919 POINTER :: box 01920 TYPE(cp_error_type), INTENT(INOUT) :: error 01921 01922 CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_gauge_new', 01923 routineP = moduleN//':'//routineN 01924 01925 INTEGER :: delta_lb(3), handle, i, iatom, ib, ibe, ibox, ibs, ie, is, j, 01926 jatom, jb, jbe, jbox, jbs, je, js, k, kb, kbe, kbox, kbs, ke, ks, 01927 lb(3), lb_local(3), natms, natms_local0, natms_local1, ng(3) 01928 REAL(dp), DIMENSION(:, :), POINTER :: r_ptr 01929 REAL(KIND=dp) :: ab, box_center(3), buf_tmp, dist, dr(3), 01930 gauge_atom_radius, offset(3), pa, pb, point(3), pra(3), r(3), res(3), 01931 summe, tmp, x, y, z 01932 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, nrm_atms_pnt 01933 REAL(kind=dp), ALLOCATABLE, 01934 DIMENSION(:, :) :: atms_pnt, ratom 01935 REAL(kind=dp), DIMENSION(:, :, :), 01936 POINTER :: grid_x, grid_y, grid_z 01937 TYPE(atomic_kind_type), DIMENSION(:), 01938 POINTER :: atomic_kind_set 01939 TYPE(cell_type), POINTER :: cell 01940 TYPE(particle_type), DIMENSION(:), 01941 POINTER :: particle_set 01942 01943 CALL timeset(routineN,handle) 01944 01945 ! 01946 CALL get_current_env(current_env=current_env,& 01947 gauge_atom_radius=gauge_atom_radius,& 01948 error=error) 01949 ! 01950 CALL get_qs_env(qs_env=qs_env,& 01951 atomic_kind_set=atomic_kind_set,& 01952 cell=cell,& 01953 particle_set=particle_set,& 01954 error=error) 01955 ! 01956 natms = SIZE(particle_set,1) 01957 dr(1) = rs_grid_x%desc%dh(1,1) 01958 dr(2) = rs_grid_x%desc%dh(2,2) 01959 dr(3) = rs_grid_x%desc%dh(3,3) 01960 lb(:) = rs_grid_x%desc%lb(:) 01961 lb_local(:) = rs_grid_x%lb_local(:) 01962 grid_x => rs_grid_x%r(:,:,:) 01963 grid_y => rs_grid_y%r(:,:,:) 01964 grid_z => rs_grid_z%r(:,:,:) 01965 ng(:) = UBOUND(grid_x) 01966 delta_lb(:) = lb_local(:)-lb(:) 01967 offset(1) = cell%hmat(1,1) * 0.5_dp 01968 offset(2) = cell%hmat(2,2) * 0.5_dp 01969 offset(3) = cell%hmat(3,3) * 0.5_dp 01970 ALLOCATE(buf(natms),ratom(3,natms),atms_pnt(3,natms),nrm_atms_pnt(natms)) 01971 ! 01972 ! find the boxes that match the grid 01973 ibs = FLOOR( REAL( delta_lb(1), dp ) * dr(1) / dbox(1) ) 01974 ibe = FLOOR( REAL( ng(1)-1+delta_lb(1), dp ) * dr(1) / dbox(1) ) 01975 jbs = FLOOR( REAL( delta_lb(2), dp ) * dr(2) / dbox(2) ) 01976 jbe = FLOOR( REAL( ng(2)-1+delta_lb(2), dp ) * dr(2) / dbox(2) ) 01977 kbs = FLOOR( REAL( delta_lb(3), dp ) * dr(3) / dbox(3) ) 01978 kbe = FLOOR( REAL( ng(3)-1+delta_lb(3), dp ) * dr(3) / dbox(3) ) 01979 ! 01980 ! go over the box-list 01981 DO kb = kbs,kbe 01982 DO jb = jbs,jbe 01983 DO ib = ibs,ibe 01984 ibox = MODULO(ib,nbox(1)) 01985 jbox = MODULO(jb,nbox(2)) 01986 kbox = MODULO(kb,nbox(3)) 01987 ! 01988 is = MAX(CEILING(REAL(ib ,dp)*dbox(1)/dr(1)), delta_lb(1))-delta_lb(1)+1 01989 ie = MIN( FLOOR(REAL(ib+1,dp)*dbox(1)/dr(1)),ng(1)-1+delta_lb(1))-delta_lb(1)+1 01990 js = MAX(CEILING(REAL(jb ,dp)*dbox(2)/dr(2)), delta_lb(2))-delta_lb(2)+1 01991 je = MIN( FLOOR(REAL(jb+1,dp)*dbox(2)/dr(2)),ng(2)-1+delta_lb(2))-delta_lb(2)+1 01992 ks = MAX(CEILING(REAL(kb ,dp)*dbox(3)/dr(3)), delta_lb(3))-delta_lb(3)+1 01993 ke = MIN( FLOOR(REAL(kb+1,dp)*dbox(3)/dr(3)),ng(3)-1+delta_lb(3))-delta_lb(3)+1 01994 ! 01995 ! sanity checks 01996 IF(.TRUE.) THEN 01997 IF(REAL(ks-1+delta_lb(3),dp)*dr(3).LT.REAL(kb,dp)*dbox(3).OR. 01998 REAL(ke-1+delta_lb(3),dp)*dr(3).GT.REAL(kb+1,dp)*dbox(3)) THEN 01999 WRITE(*,*) 'box_k',REAL(kb,dp)*dbox(3),REAL(kb+1,dp)*dbox(3) 02000 WRITE(*,*) 'point_k',REAL(ks-1+delta_lb(3),dp)*dr(3),REAL(ke-1+delta_lb(3),dp)*dr(3) 02001 WRITE(*,*) 'ibox',ibox,'jbox',jbox,'kbox',kbox 02002 WRITE(*,*) 'is,ie',is,ie,' js,je',js,je,' ks,ke',ks,ke 02003 WRITE(*,*) 'ibs,ibe',ibs,ibe,' jbs,jbe',jbs,jbe,' kbs,kbe',kbs,kbe 02004 STOP 'we stop_k' 02005 ENDIF 02006 IF(REAL(js-1+delta_lb(2),dp)*dr(2).LT.REAL(jb,dp)*dbox(2).OR. 02007 REAL(je-1+delta_lb(2),dp)*dr(2).GT.REAL(jb+1,dp)*dbox(2)) THEN 02008 WRITE(*,*) 'box_j',REAL(jb,dp)*dbox(2),REAL(jb+1,dp)*dbox(2) 02009 WRITE(*,*) 'point_j',REAL(js-1+delta_lb(2),dp)*dr(2),REAL(je-1+delta_lb(2),dp)*dr(2) 02010 WRITE(*,*) 'is,ie',is,ie,' js,je',js,je,' ks,ke',ks,ke 02011 WRITE(*,*) 'ibs,ibe',ibs,ibe,' jbs,jbe',jbs,jbe,' kbs,kbe',kbs,kbe 02012 STOP 'we stop_j' 02013 ENDIF 02014 IF(REAL(is-1+delta_lb(1),dp)*dr(1).LT.REAL(ib,dp)*dbox(1).OR. 02015 REAL(ie-1+delta_lb(1),dp)*dr(1).GT.REAL(ib+1,dp)*dbox(1)) THEN 02016 WRITE(*,*) 'box_i',REAL(ib,dp)*dbox(1),REAL(ib+1,dp)*dbox(1) 02017 WRITE(*,*) 'point_i',REAL(is-1+delta_lb(1),dp)*dr(1),REAL(ie-1+delta_lb(1),dp)*dr(1) 02018 WRITE(*,*) 'is,ie',is,ie,' js,je',js,je,' ks,ke',ks,ke 02019 WRITE(*,*) 'ibs,ibe',ibs,ibe,' jbs,jbe',jbs,jbe,' kbs,kbe',kbs,kbe 02020 STOP 'we stop_i' 02021 ENDIF 02022 ENDIF 02023 ! 02024 ! the center of the box 02025 box_center(1) = ( REAL(ibox,dp) + 0.5_dp ) * dbox(1) 02026 box_center(2) = ( REAL(jbox,dp) + 0.5_dp ) * dbox(2) 02027 box_center(3) = ( REAL(kbox,dp) + 0.5_dp ) * dbox(3) 02028 ! 02029 ! find the atoms that are in the overlaping boxes 02030 natms_local0 = box(ibox,jbox,kbox)%n 02031 r_ptr => box(ibox,jbox,kbox)%r 02032 ! 02033 ! go over the grid inside the box 02034 IF(natms_local0.GT.0) THEN 02035 ! 02036 ! here there are some atoms... 02037 DO k = ks,ke 02038 DO j = js,je 02039 DO i = is,ie 02040 point(1) = REAL(i-1+delta_lb(1),dp)*dr(1) 02041 point(2) = REAL(j-1+delta_lb(2),dp)*dr(2) 02042 point(3) = REAL(k-1+delta_lb(3),dp)*dr(3) 02043 point = pbc(point,cell) 02044 ! 02045 ! compute atom-point distances 02046 natms_local1 = 0 02047 DO iatom = 1,natms_local0 02048 r(:) = pbc(r_ptr(:,iatom)-point(:),cell) + point(:)!needed? 02049 dist = (r(1)-point(1))**2 + (r(2)-point(2))**2 + (r(3)-point(3))**2 02050 IF(dist.LT.gauge_atom_radius**2) THEN 02051 natms_local1 = natms_local1 + 1 02052 ratom(:,natms_local1) = r(:) 02053 ! 02054 ! compute the distance atoms-point 02055 x = point(1) - r(1) 02056 y = point(2) - r(2) 02057 z = point(3) - r(3) 02058 atms_pnt(1,natms_local1) = x 02059 atms_pnt(2,natms_local1) = y 02060 atms_pnt(3,natms_local1) = z 02061 nrm_atms_pnt(natms_local1) = SQRT( x*x + y*y + z*z ) 02062 ENDIF 02063 ENDDO 02064 ! 02065 ! 02066 IF(natms_local1.GT.0) THEN 02067 ! 02068 ! build the step 02069 DO iatom = 1,natms_local1 02070 buf_tmp = 1.0_dp 02071 pra(1) = atms_pnt(1,iatom) 02072 pra(2) = atms_pnt(2,iatom) 02073 pra(3) = atms_pnt(3,iatom) 02074 pa = nrm_atms_pnt(iatom) 02075 DO jatom = 1,natms_local1 02076 IF(iatom.EQ.jatom) CYCLE 02077 pb = nrm_atms_pnt(jatom) 02078 x = pra(1) - atms_pnt(1,jatom) 02079 y = pra(2) - atms_pnt(2,jatom) 02080 z = pra(3) - atms_pnt(3,jatom) 02081 ab = SQRT( x*x + y*y + z*z ) 02082 ! 02083 tmp = ( pa - pb ) / ab 02084 tmp = 0.5_dp * ( 3.0_dp - tmp * tmp ) * tmp 02085 tmp = 0.5_dp * ( 3.0_dp - tmp * tmp ) * tmp 02086 tmp = 0.5_dp * ( 3.0_dp - tmp * tmp ) * tmp 02087 buf_tmp = buf_tmp * 0.5_dp * ( 1.0_dp - tmp ) 02088 ENDDO 02089 buf(iatom) = buf_tmp 02090 ENDDO 02091 res(1) = 0.0_dp 02092 res(2) = 0.0_dp 02093 res(3) = 0.0_dp 02094 summe = 0.0_dp 02095 DO iatom = 1,natms_local1 02096 res(1) = res(1) + ratom(1,iatom) * buf(iatom) 02097 res(2) = res(2) + ratom(2,iatom) * buf(iatom) 02098 res(3) = res(3) + ratom(3,iatom) * buf(iatom) 02099 summe = summe + buf(iatom) 02100 ENDDO 02101 res(1) = res(1) / summe 02102 res(2) = res(2) / summe 02103 res(3) = res(3) / summe 02104 grid_x(i,j,k) = point(1) - res(1) 02105 grid_y(i,j,k) = point(2) - res(2) 02106 grid_z(i,j,k) = point(3) - res(3) 02107 ELSE 02108 grid_x(i,j,k) = 0.0_dp 02109 grid_y(i,j,k) = 0.0_dp 02110 grid_z(i,j,k) = 0.0_dp 02111 ENDIF 02112 ENDDO ! grid 02113 ENDDO 02114 ENDDO 02115 ! 02116 ELSE 02117 ! 02118 ! here there is no atom 02119 DO k = ks,ke 02120 DO j = js,je 02121 DO i = is,ie 02122 grid_x(i,j,k) = 0.0_dp 02123 grid_y(i,j,k) = 0.0_dp 02124 grid_z(i,j,k) = 0.0_dp 02125 ENDDO ! grid 02126 ENDDO 02127 ENDDO 02128 ! 02129 ENDIF 02130 ! 02131 ENDDO ! list 02132 ENDDO 02133 ENDDO 02134 02135 DEALLOCATE(buf,ratom,atms_pnt,nrm_atms_pnt) 02136 02137 CALL timestop(handle) 02138 02139 END SUBROUTINE collocate_gauge_new 02140 02141 SUBROUTINE deallocate_box(box) 02142 TYPE(box_type), DIMENSION(:, :, :), 02143 POINTER :: box 02144 02145 INTEGER :: i, j, k 02146 02147 IF(ASSOCIATED(box)) THEN 02148 DO k = LBOUND(box,3),UBOUND(box,3) 02149 DO j = LBOUND(box,2),UBOUND(box,2) 02150 DO i = LBOUND(box,1),UBOUND(box,1) 02151 IF(ASSOCIATED(box(i,j,k)%r)) THEN 02152 DEALLOCATE(box(i,j,k)%r) 02153 ENDIF 02154 ENDDO 02155 ENDDO 02156 ENDDO 02157 DEALLOCATE(box) 02158 ENDIF 02159 END SUBROUTINE deallocate_box 02160 END SUBROUTINE current_set_gauge 02161 02162 02163 SUBROUTINE current_build_chi(current_env,qs_env,iB,error) 02164 ! 02165 TYPE(current_env_type) :: current_env 02166 TYPE(qs_environment_type), POINTER :: qs_env 02167 INTEGER, INTENT(IN) :: iB 02168 TYPE(cp_error_type), INTENT(INOUT) :: error 02169 02170 IF(current_env%full) THEN 02171 CALL current_build_chi_many_centers(current_env,qs_env,iB,error) 02172 ELSE 02173 CALL current_build_chi_one_center(current_env,qs_env,iB,error) 02174 ENDIF 02175 02176 END SUBROUTINE current_build_chi 02177 02178 02179 SUBROUTINE current_build_chi_many_centers(current_env,qs_env,iB,error) 02180 ! 02181 TYPE(current_env_type) :: current_env 02182 TYPE(qs_environment_type), POINTER :: qs_env 02183 INTEGER, INTENT(IN) :: iB 02184 TYPE(cp_error_type), INTENT(INOUT) :: error 02185 02186 CHARACTER(LEN=*), PARAMETER :: 02187 routineN = 'current_build_chi_many_centers', 02188 routineP = moduleN//':'//routineN 02189 02190 INTEGER :: handle, icenter, idir, idir2, ii, iiB, iii, iiiB, ispin, 02191 istat, istate, j, jstate, max_states, nao, natom, nbr_center(2), nmo, 02192 nspins, nstate_loc, nstates(2), output_unit 02193 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 02194 INTEGER, DIMENSION(:), POINTER :: rbs 02195 LOGICAL :: chi_pbc, failure, gapw 02196 REAL(dp) :: chi(3), chi_tmp, contrib, 02197 contrib2, dk(3), 02198 int_current(3), 02199 int_current_tmp, maxocc 02200 TYPE(array_i1d_obj) :: row_blk_sizes 02201 TYPE(cell_type), POINTER :: cell 02202 TYPE(cp_2d_i_p_type), DIMENSION(:), 02203 POINTER :: center_list 02204 TYPE(cp_2d_r_p_type), DIMENSION(:), 02205 POINTER :: centers_set 02206 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02207 POINTER :: op_mom_ao, op_p_ao 02208 TYPE(cp_dbcsr_p_type), DIMENSION(:, :), 02209 POINTER :: op_mom_der_ao 02210 TYPE(cp_fm_p_type), DIMENSION(:), 02211 POINTER :: p_rxp, psi0_order, r_p1, r_p2 02212 TYPE(cp_fm_p_type), DIMENSION(:, :), 02213 POINTER :: psi1_D, psi1_p, psi1_rxp, 02214 rr_p1, rr_p2, rr_rxp 02215 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 02216 TYPE(cp_fm_type), POINTER :: mo_coeff, psi0, psi_D, 02217 psi_p1, psi_p2, psi_rxp 02218 TYPE(cp_logger_type), POINTER :: logger 02219 TYPE(cp_para_env_type), POINTER :: para_env 02220 TYPE(dbcsr_distribution_obj), POINTER :: dbcsr_dist 02221 TYPE(dft_control_type), POINTER :: dft_control 02222 TYPE(mo_set_p_type), DIMENSION(:), 02223 POINTER :: mos 02224 TYPE(neighbor_list_set_p_type), 02225 DIMENSION(:), POINTER :: sab_all, sab_orb 02226 TYPE(particle_type), DIMENSION(:), 02227 POINTER :: particle_set 02228 02229 failure = .FALSE. 02230 ! 02231 CALL timeset(routineN,handle) 02232 ! 02233 NULLIFY(dft_control,mos,para_env,mo_coeff,op_mom_ao,& 02234 & op_mom_der_ao,center_list,centers_set,psi0,psi_rxp,psi_D,psi_p1,psi_p2,& 02235 & op_p_ao,psi1_p,psi1_rxp,psi1_D,p_rxp,r_p1,r_p2,rr_rxp,rr_p1,rr_p2,& 02236 & cell,psi0_order,particle_set) 02237 02238 logger => cp_error_get_logger(error) 02239 output_unit= cp_logger_get_default_io_unit(logger) 02240 02241 CALL get_qs_env(qs_env=qs_env,& 02242 & dft_control=dft_control,& 02243 & mos=mos,& 02244 & para_env=para_env,& 02245 & cell=cell,& 02246 & dbcsr_dist=dbcsr_dist,& 02247 & particle_set=particle_set,& 02248 & sab_all=sab_all,& 02249 & sab_orb=sab_orb,& 02250 & error=error) 02251 02252 nspins = dft_control%nspins 02253 gapw = dft_control%qs_control%gapw 02254 02255 CALL get_current_env(current_env=current_env,& 02256 & chi_pbc=chi_pbc,& 02257 & nao=nao,& 02258 & nbr_center=nbr_center,& 02259 & center_list=center_list,& 02260 & centers_set=centers_set,& 02261 & psi1_p=psi1_p,& 02262 & psi1_rxp=psi1_rxp,& 02263 & psi1_D=psi1_D,& 02264 & nstates=nstates,& 02265 & psi0_order=psi0_order,& 02266 & error=error) 02267 ! 02268 ! get max nbr of states per center 02269 max_states = 0 02270 DO ispin = 1,nspins 02271 DO icenter = 1,nbr_center(ispin) 02272 max_states = MAX(max_states,center_list(ispin)%array(1,icenter+1)& 02273 & -center_list(ispin)%array(1,icenter)) 02274 ENDDO 02275 ENDDO 02276 ! 02277 ! Allocate sparse matrices for dipole, quadrupole and their derivatives => 9x3 02278 ! Remember the derivatives are antisymmetric 02279 CALL cp_dbcsr_allocate_matrix_set(op_mom_ao,9,error=error) 02280 CALL cp_dbcsr_allocate_matrix_set(op_mom_der_ao,9,3,error=error) 02281 ! 02282 ! prepare for allocation 02283 natom = SIZE(particle_set,1) 02284 ALLOCATE (first_sgf(natom),STAT=istat) 02285 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02286 ALLOCATE (last_sgf(natom),STAT=istat) 02287 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02288 CALL get_particle_set(particle_set=particle_set,& 02289 first_sgf=first_sgf,& 02290 last_sgf=last_sgf,error=error) 02291 ALLOCATE (rbs(natom), STAT=istat) 02292 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02293 CALL convert_offsets_to_sizes (first_sgf, rbs, last_sgf) 02294 CALL array_nullify (row_blk_sizes) 02295 CALL array_new (row_blk_sizes, rbs, gift=.TRUE.) 02296 DEALLOCATE (first_sgf,STAT=istat) 02297 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02298 DEALLOCATE (last_sgf,STAT=istat) 02299 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02300 ! 02301 ! 02302 ALLOCATE(op_mom_ao(1)%matrix) 02303 CALL cp_dbcsr_init(op_mom_ao(1)%matrix,error=error) 02304 CALL cp_dbcsr_create(matrix=op_mom_ao(1)%matrix, & 02305 name="op_mom", & 02306 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry,& 02307 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 02308 nblks=0, nze=0, mutable_work=.TRUE., & 02309 error=error) 02310 CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix,sab_all,error=error) 02311 02312 DO idir2=1,3 02313 ALLOCATE(op_mom_der_ao(1,idir2)%matrix) 02314 CALL cp_dbcsr_init(op_mom_der_ao(1,idir2)%matrix, error=error) 02315 CALL cp_dbcsr_copy(op_mom_der_ao(1,idir2)%matrix,op_mom_ao(1)%matrix,& 02316 "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir2))),error) 02317 ENDDO 02318 02319 DO idir = 2,SIZE(op_mom_ao,1) 02320 ALLOCATE(op_mom_ao(idir)%matrix) 02321 CALL cp_dbcsr_init(op_mom_ao(idir)%matrix, error=error) 02322 CALL cp_dbcsr_copy(op_mom_ao(idir)%matrix,op_mom_ao(1)%matrix,& 02323 "op_mom_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))),error) 02324 DO idir2=1,3 02325 ALLOCATE(op_mom_der_ao(idir,idir2)%matrix) 02326 CALL cp_dbcsr_init(op_mom_der_ao(idir,idir2)%matrix, error=error) 02327 CALL cp_dbcsr_copy(op_mom_der_ao(idir,idir2)%matrix,op_mom_ao(1)%matrix,& 02328 "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir*idir2))),error) 02329 ENDDO 02330 ENDDO 02331 ! 02332 CALL cp_dbcsr_allocate_matrix_set(op_p_ao,3,error=error) 02333 ALLOCATE(op_p_ao(1)%matrix) 02334 CALL cp_dbcsr_init(op_p_ao(1)%matrix,error=error) 02335 CALL cp_dbcsr_create(matrix=op_p_ao(1)%matrix, & 02336 name="op_p_ao", & 02337 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric,& 02338 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 02339 nblks=0, nze=0, mutable_work=.TRUE., & 02340 error=error) 02341 CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb, error=error) 02342 02343 DO idir=2,3 02344 ALLOCATE(op_p_ao(idir)%matrix) 02345 CALL cp_dbcsr_init(op_p_ao(idir)%matrix, error=error) 02346 CALL cp_dbcsr_copy(op_p_ao(idir)%matrix,op_p_ao(1)%matrix,& 02347 "op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))),error) 02348 ENDDO 02349 ! 02350 ! 02351 CALL array_release (row_blk_sizes) 02352 ! 02353 ! 02354 ! Allocate full matrices for only one vector 02355 mo_coeff => psi0_order(1)%matrix 02356 NULLIFY(tmp_fm_struct) 02357 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,& 02358 & ncol_global=max_states,para_env=para_env,& 02359 & context=mo_coeff%matrix_struct%context,& 02360 & error=error) 02361 CALL cp_fm_create(psi0,tmp_fm_struct,error=error) 02362 CALL cp_fm_create(psi_D,tmp_fm_struct,error=error) 02363 CALL cp_fm_create(psi_rxp,tmp_fm_struct,error=error) 02364 CALL cp_fm_create(psi_p1,tmp_fm_struct,error=error) 02365 CALL cp_fm_create(psi_p2,tmp_fm_struct,error=error) 02366 CALL cp_fm_struct_release(tmp_fm_struct,error=error) 02367 ! 02368 ALLOCATE(p_rxp(3),r_p1(3),r_p2(3),rr_rxp(9,3),rr_p1(9,3),rr_p2(9,3),STAT=istat) 02369 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02370 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,& 02371 & ncol_global=max_states,para_env=para_env,& 02372 & context=mo_coeff%matrix_struct%context,& 02373 & error=error) 02374 DO idir = 1,3 02375 CALL cp_fm_create(p_rxp(idir)%matrix,tmp_fm_struct,error=error) 02376 CALL cp_fm_create(r_p1(idir)%matrix,tmp_fm_struct,error=error) 02377 CALL cp_fm_create(r_p2(idir)%matrix,tmp_fm_struct,error=error) 02378 DO idir2 = 1,9 02379 CALL cp_fm_create(rr_rxp(idir2,idir)%matrix,tmp_fm_struct,error=error) 02380 CALL cp_fm_create(rr_p1(idir2,idir)%matrix,tmp_fm_struct,error=error) 02381 CALL cp_fm_create(rr_p2(idir2,idir)%matrix,tmp_fm_struct,error=error) 02382 ENDDO 02383 ENDDO 02384 CALL cp_fm_struct_release(tmp_fm_struct,error=error) 02385 ! 02386 ! 02387 ! 02388 ! recompute the linear momentum matrices 02389 CALL build_lin_mom_matrix(qs_env,op_p_ao,error) 02390 !CALL p_xyz_ao(op_p_ao,qs_env,minimum_image=.FALSE.,error=error) 02391 ! 02392 ! 02393 ! get iiB and iiiB 02394 CALL set_vecp(iB,iiB,iiiB) 02395 DO ispin = 1,nspins 02396 ! 02397 ! get ground state MOS 02398 nmo = nstates(ispin) 02399 mo_coeff => psi0_order(ispin)%matrix 02400 CALL get_mo_set(mo_set=mos(ispin)%mo_set,maxocc=maxocc) 02401 ! 02402 ! Initialize the temporary vector chi 02403 chi = 0.0_dp 02404 int_current = 0.0_dp 02405 ! 02406 ! Start loop over the occupied states 02407 DO icenter = 1,nbr_center(ispin) 02408 ! 02409 ! Get the Wannier center of the istate-th ground state orbital 02410 dk(1:3) = centers_set(ispin)%array(1:3,icenter) 02411 ! 02412 ! Compute the multipole integrals for the state istate, 02413 ! using as reference center the corresponding Wannier center 02414 DO idir = 1,9 02415 CALL cp_dbcsr_set(op_mom_ao(idir)%matrix,0.0_dp,error=error) 02416 DO idir2 = 1,3 02417 CALL cp_dbcsr_set(op_mom_der_ao(idir,idir2)%matrix,0.0_dp,error=error) 02418 ENDDO 02419 ENDDO 02420 CALL rRc_xyz_der_ao(op_mom_ao,op_mom_der_ao,qs_env,dk,order=2,& 02421 & minimum_image=.FALSE.,soft=gapw,error=error) 02422 ! 02423 ! collecte the states that belong to a given center 02424 CALL cp_fm_set_all(psi0,0.0_dp,error=error) 02425 CALL cp_fm_set_all(psi_rxp,0.0_dp,error=error) 02426 CALL cp_fm_set_all(psi_D,0.0_dp,error=error) 02427 CALL cp_fm_set_all(psi_p1,0.0_dp,error=error) 02428 CALL cp_fm_set_all(psi_p2,0.0_dp,error=error) 02429 nstate_loc = center_list(ispin)%array(1,icenter+1)-center_list(ispin)%array(1,icenter) 02430 jstate = 1 02431 DO j = center_list(ispin)%array(1,icenter),center_list(ispin)%array(1,icenter+1)-1 02432 istate = center_list(ispin)%array(2,j) 02433 ! 02434 ! block the states that belong to this center 02435 CALL cp_fm_to_fm(mo_coeff,psi0,1,istate,jstate) 02436 ! 02437 CALL cp_fm_to_fm(psi1_rxp(ispin,iB)%matrix,psi_rxp,1,istate,jstate) 02438 IF(current_env%full) CALL cp_fm_to_fm(psi1_D(ispin,iB)%matrix,psi_D,1,istate,jstate) 02439 ! 02440 ! psi1_p_iiB_istate and psi1_p_iiiB_istate 02441 CALL cp_fm_to_fm(psi1_p(ispin, iiB)%matrix,psi_p1,1,istate,jstate) 02442 CALL cp_fm_to_fm(psi1_p(ispin,iiiB)%matrix,psi_p2,1,istate,jstate) 02443 ! 02444 jstate = jstate+1 02445 ENDDO ! istate 02446 ! 02447 ! scale the ordered mos 02448 IF(current_env%full) CALL cp_fm_scale_and_add(1.0_dp,psi_rxp,-1.0_dp,psi_D,error=error) 02449 ! 02450 DO idir = 1,3 02451 CALL set_vecp(idir,ii,iii) 02452 CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix,psi_rxp,& 02453 & p_rxp(idir)%matrix,ncol=nstate_loc,alpha=1.e0_dp,& 02454 & error=error) 02455 IF(iiiB.EQ.iii.OR.iiiB.EQ.ii) THEN 02456 CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix,psi_p1,& 02457 & r_p1(idir)%matrix,ncol=nstate_loc,alpha=1.e0_dp,& 02458 & error=error) 02459 ENDIF 02460 IF(iiB.EQ.iii.OR.iiB.EQ.ii) THEN 02461 CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix,psi_p2,& 02462 & r_p2(idir)%matrix,ncol=nstate_loc,alpha=1.e0_dp,& 02463 & error=error) 02464 ENDIF 02465 DO idir2 = 1,9 02466 IF(idir2.EQ.ii.OR.idir2.EQ.iii) THEN 02467 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2,idir)%matrix,psi_rxp,& 02468 & rr_rxp(idir2,idir)%matrix,ncol=nstate_loc,alpha=1.e0_dp,& 02469 & error=error) 02470 ENDIF 02471 ! 02472 IF(idir2.EQ.ind_m2(ii,iiiB).OR.idir2.EQ.ind_m2(iii,iiiB).OR.idir2.EQ.iiiB) THEN 02473 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2,idir)%matrix,psi_p1,& 02474 & rr_p1(idir2,idir)%matrix,ncol=nstate_loc,alpha=1.e0_dp,& 02475 & error=error) 02476 ENDIF 02477 ! 02478 IF(idir2.EQ.ind_m2(ii,iiB).OR.idir2.EQ.ind_m2(iii,iiB).OR.idir2.EQ.iiB) THEN 02479 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2,idir)%matrix,psi_p2,& 02480 & rr_p2(idir2,idir)%matrix,ncol=nstate_loc,alpha=1.e0_dp,& 02481 & error=error) 02482 ENDIF 02483 ENDDO 02484 ENDDO 02485 ! 02486 ! Multuply left and right by the appropriate coefficients and sum into the 02487 ! correct component of the chi tensor using the appropriate multiplicative factor 02488 ! (don't forget the occupation number) 02489 ! Loop over the cartesian components of the tensor 02490 ! The loop over the components of the external field is external, thereby 02491 ! only one column of the chi tensor is computed here 02492 DO idir = 1,3 02493 chi_tmp = 0.0_dp 02494 int_current_tmp = 0.0_dp 02495 ! 02496 ! get ii and iii 02497 CALL set_vecp(idir,ii,iii) 02498 ! 02499 ! term: 2[C0| (r-dk)_ii |d_iii(C1(rxp-D))]-2[C0| (r-dk)_iii |d_ii(C1(rxp-D))] 02500 ! the factor 2 should be already included in the matrix elements 02501 contrib = 0.0_dp 02502 CALL cp_fm_trace(psi0,rr_rxp(ii,iii)%matrix,contrib,error=error) 02503 chi_tmp = chi_tmp + 2.0_dp * contrib 02504 ! 02505 contrib = 0.0_dp 02506 CALL cp_fm_trace(psi0,rr_rxp(iii,ii)%matrix,contrib,error=error) 02507 chi_tmp = chi_tmp - 2.0_dp * contrib 02508 ! 02509 ! correction: dk_ii*2[C0| d_iii(C1(rxp-D))] - dk_iii*2[C0| d_ii(C1(rxp-D))] 02510 ! factor 2 not included in the matrix elements 02511 contrib = 0.0_dp 02512 CALL cp_fm_trace(psi0,p_rxp(iii)%matrix,contrib,error=error) 02513 IF(.NOT.chi_pbc) chi_tmp = chi_tmp + 2.0_dp * dk(ii) * contrib 02514 int_current_tmp = int_current_tmp + 2.0_dp * contrib 02515 ! 02516 contrib2 = 0.0_dp 02517 CALL cp_fm_trace(psi0,p_rxp(ii)%matrix,contrib2,error=error) 02518 IF(.NOT.chi_pbc) chi_tmp = chi_tmp - 2.0_dp * dk(iii) * contrib2 02519 ! 02520 ! term: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] \ 02521 ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))] 02522 ! the factor 2 should be already included in the matrix elements 02523 contrib = 0.0_dp 02524 idir2 = ind_m2(ii,iiB) 02525 CALL cp_fm_trace(psi0,rr_p2(idir2,iii)%matrix,contrib,error=error) 02526 chi_tmp = chi_tmp - 2.0_dp * contrib 02527 contrib2 = 0.0_dp 02528 IF(iiB==iii) THEN 02529 CALL cp_fm_trace(psi0,r_p2(ii)%matrix,contrib2,error=error) 02530 chi_tmp = chi_tmp - contrib2 02531 ENDIF 02532 ! 02533 contrib = 0.0_dp 02534 idir2 = ind_m2(iii,iiB) 02535 CALL cp_fm_trace(psi0,rr_p2(idir2,ii)%matrix,contrib,error=error) 02536 chi_tmp = chi_tmp + 2.0_dp * contrib 02537 contrib2 = 0.0_dp 02538 IF(iiB==ii) THEN 02539 CALL cp_fm_trace(psi0,r_p2(iii)%matrix,contrib2,error=error) 02540 chi_tmp = chi_tmp + contrib2 02541 ENDIF 02542 ! 02543 ! correction: -dk_ii * 2[C0|(r-dk)_iiB | d_iii(C1(piiiB))] \ 02544 ! +dk_iii * 2[C0|(r-dk)_iiB | d_ii(C1(piiiB))] 02545 ! the factor 2 should be already included in the matrix elements 02546 ! no additional correction terms because of the orthogonality between C0 and C1 02547 contrib = 0.0_dp 02548 CALL cp_fm_trace(psi0,rr_p2(iiB,iii)%matrix,contrib,error=error) 02549 IF(.NOT.chi_pbc) chi_tmp = chi_tmp - 2.0_dp * dk(ii) * contrib 02550 int_current_tmp = int_current_tmp - 2.0_dp * contrib 02551 ! 02552 contrib2 = 0.0_dp 02553 CALL cp_fm_trace(psi0,rr_p2(iiB,ii)%matrix,contrib2,error=error) 02554 IF(.NOT.chi_pbc) chi_tmp = chi_tmp + 2.0_dp * dk(iii) * contrib2 02555 ! 02556 ! term: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] \ 02557 ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))] 02558 ! the factor 2 should be already included in the matrix elements 02559 contrib = 0.0_dp 02560 idir2 = ind_m2(ii,iiiB) 02561 CALL cp_fm_trace(psi0,rr_p1(idir2,iii)%matrix,contrib,error=error) 02562 chi_tmp = chi_tmp + 2.0_dp * contrib 02563 contrib2 = 0.0_dp 02564 IF(iiiB==iii) THEN 02565 CALL cp_fm_trace(psi0,r_p1(ii)%matrix,contrib2,error=error) 02566 chi_tmp = chi_tmp + contrib2 02567 ENDIF 02568 ! 02569 contrib = 0.0_dp 02570 idir2 = ind_m2(iii,iiiB) 02571 CALL cp_fm_trace(psi0,rr_p1(idir2,ii)%matrix,contrib,error=error) 02572 chi_tmp = chi_tmp - 2.0_dp * contrib 02573 contrib2 = 0.0_dp 02574 IF(iiiB==ii) THEN 02575 CALL cp_fm_trace(psi0,r_p1(iii)%matrix,contrib2,error=error) 02576 chi_tmp = chi_tmp - contrib2 02577 ENDIF 02578 ! 02579 ! correction: +dk_ii * 2[C0|(r-dk)_iiiB | d_iii(C1(piiB))] +\ 02580 ! -dk_iii * 2[C0|(r-dk)_iiiB | d_ii(C1(piiB))] 02581 ! the factor 2 should be already included in the matrix elements 02582 contrib = 0.0_dp 02583 CALL cp_fm_trace(psi0,rr_p1(iiiB,iii)%matrix,contrib,error=error) 02584 IF(.NOT.chi_pbc) chi_tmp = chi_tmp + 2.0_dp * dk(ii) * contrib 02585 int_current_tmp = int_current_tmp + 2.0_dp * contrib 02586 ! 02587 contrib2 = 0.0_dp 02588 CALL cp_fm_trace(psi0,rr_p1(iiiB,ii)%matrix,contrib2,error=error) 02589 IF(.NOT.chi_pbc) chi_tmp = chi_tmp - 2.0_dp * dk(iii) * contrib2 02590 ! 02591 ! accumulate 02592 chi(idir) = chi(idir) + maxocc * chi_tmp 02593 int_current(iii) = int_current(iii) + int_current_tmp 02594 ENDDO ! idir 02595 02596 ENDDO ! icenter 02597 ! 02598 DO idir = 1,3 02599 current_env%chi_tensor(idir,iB,ispin) = current_env%chi_tensor(idir,iB,ispin)+& 02600 & chi(idir) 02601 IF (output_unit>0) THEN 02602 !WRITE(output_unit,'(A,E12.6)') ' chi_'//ACHAR(119+idir)//ACHAR(119+iB)//& 02603 ! & ' = ',chi(idir) 02604 !WRITE(output_unit,'(A,E12.6)') ' analytic \int j_'//ACHAR(119+idir)//ACHAR(119+iB)//& 02605 ! & '(r) d^3r = ',int_current(idir) 02606 ENDIF 02607 ENDDO 02608 ! 02609 ENDDO ! ispin 02610 ! 02611 ! deallocate the sparse matrices 02612 CALL cp_dbcsr_deallocate_matrix_set(op_mom_ao,error=error) 02613 CALL cp_dbcsr_deallocate_matrix_set(op_mom_der_ao,error=error) 02614 CALL cp_dbcsr_deallocate_matrix_set(op_p_ao,error=error) 02615 02616 CALL cp_fm_release(psi0,error=error) 02617 CALL cp_fm_release(psi_rxp,error=error) 02618 CALL cp_fm_release(psi_D,error=error) 02619 CALL cp_fm_release(psi_p1,error=error) 02620 CALL cp_fm_release(psi_p2,error=error) 02621 DO idir = 1,3 02622 CALL cp_fm_release(p_rxp(idir)%matrix,error=error) 02623 CALL cp_fm_release(r_p1(idir)%matrix,error=error) 02624 CALL cp_fm_release(r_p2(idir)%matrix,error=error) 02625 DO idir2 = 1,9 02626 CALL cp_fm_release(rr_rxp(idir2,idir)%matrix,error=error) 02627 CALL cp_fm_release(rr_p1(idir2,idir)%matrix,error=error) 02628 CALL cp_fm_release(rr_p2(idir2,idir)%matrix,error=error) 02629 ENDDO 02630 ENDDO 02631 DEALLOCATE(p_rxp,r_p1,r_p2,rr_rxp,rr_p1,rr_p2,STAT=istat) 02632 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02633 02634 CALL timestop(handle) 02635 02636 END SUBROUTINE current_build_chi_many_centers 02637 02638 02639 SUBROUTINE current_build_chi_one_center(current_env,qs_env,iB,error) 02640 ! 02641 TYPE(current_env_type) :: current_env 02642 TYPE(qs_environment_type), POINTER :: qs_env 02643 INTEGER, INTENT(IN) :: iB 02644 TYPE(cp_error_type), INTENT(INOUT) :: error 02645 02646 CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_chi_one_center', 02647 routineP = moduleN//':'//routineN 02648 02649 INTEGER :: handle, idir, idir2, iiB, iiiB, ispin, istat, jdir, jjdir, 02650 kdir, max_states, nao, natom, nbr_center(2), nmo, nspins, nstates(2), 02651 output_unit 02652 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 02653 INTEGER, DIMENSION(:), POINTER :: rbs 02654 LOGICAL :: chi_pbc, failure, gapw 02655 REAL(dp) :: chi(3), contrib, dk(3), 02656 int_current(3), maxocc 02657 TYPE(array_i1d_obj) :: row_blk_sizes 02658 TYPE(cell_type), POINTER :: cell 02659 TYPE(cp_2d_i_p_type), DIMENSION(:), 02660 POINTER :: center_list 02661 TYPE(cp_2d_r_p_type), DIMENSION(:), 02662 POINTER :: centers_set 02663 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02664 POINTER :: op_mom_ao, op_p_ao 02665 TYPE(cp_dbcsr_p_type), DIMENSION(:, :), 02666 POINTER :: op_mom_der_ao 02667 TYPE(cp_fm_p_type), DIMENSION(:), 02668 POINTER :: psi0_order 02669 TYPE(cp_fm_p_type), DIMENSION(:, :), 02670 POINTER :: psi1_D, psi1_p, psi1_rxp 02671 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 02672 TYPE(cp_fm_type), POINTER :: buf, mo_coeff 02673 TYPE(cp_logger_type), POINTER :: logger 02674 TYPE(cp_para_env_type), POINTER :: para_env 02675 TYPE(dbcsr_distribution_obj), POINTER :: dbcsr_dist 02676 TYPE(dft_control_type), POINTER :: dft_control 02677 TYPE(mo_set_p_type), DIMENSION(:), 02678 POINTER :: mos 02679 TYPE(neighbor_list_set_p_type), 02680 DIMENSION(:), POINTER :: sab_all, sab_orb 02681 TYPE(particle_type), DIMENSION(:), 02682 POINTER :: particle_set 02683 02684 failure = .FALSE. 02685 ! 02686 CALL timeset(routineN,handle) 02687 ! 02688 NULLIFY(dft_control,mos,para_env,mo_coeff,op_mom_ao,& 02689 & op_mom_der_ao,center_list,centers_set,& 02690 & op_p_ao,psi1_p,psi1_rxp,psi1_D,buf,cell,psi0_order) 02691 02692 logger => cp_error_get_logger(error) 02693 output_unit= cp_logger_get_default_io_unit(logger) 02694 02695 CALL get_qs_env(qs_env=qs_env,& 02696 & dft_control=dft_control,& 02697 & mos=mos,& 02698 & para_env=para_env,& 02699 & cell=cell,& 02700 & particle_set=particle_set,& 02701 & sab_all=sab_all,& 02702 & sab_orb=sab_orb,& 02703 & dbcsr_dist=dbcsr_dist,& 02704 & error=error) 02705 02706 nspins = dft_control%nspins 02707 gapw = dft_control%qs_control%gapw 02708 02709 CALL get_current_env(current_env=current_env,& 02710 & chi_pbc=chi_pbc,& 02711 & nao=nao,& 02712 & nbr_center=nbr_center,& 02713 & center_list=center_list,& 02714 & centers_set=centers_set,& 02715 & psi1_p=psi1_p,& 02716 & psi1_rxp=psi1_rxp,& 02717 & psi1_D=psi1_D,& 02718 & nstates=nstates,& 02719 & psi0_order=psi0_order,& 02720 & error=error) 02721 ! 02722 max_states = MAXVAL(nstates(1:nspins)) 02723 ! 02724 ! Allocate sparse matrices for dipole, quadrupole and their derivatives => 9x3 02725 ! Remember the derivatives are antisymmetric 02726 CALL cp_dbcsr_allocate_matrix_set(op_mom_ao,9,error=error) 02727 CALL cp_dbcsr_allocate_matrix_set(op_mom_der_ao,9,3,error=error) 02728 ! 02729 ! prepare for allocation 02730 natom = SIZE(particle_set,1) 02731 ALLOCATE (first_sgf(natom),STAT=istat) 02732 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02733 ALLOCATE (last_sgf(natom),STAT=istat) 02734 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02735 CALL get_particle_set(particle_set=particle_set,& 02736 first_sgf=first_sgf,& 02737 last_sgf=last_sgf,error=error) 02738 ALLOCATE (rbs(natom), STAT=istat) 02739 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02740 CALL convert_offsets_to_sizes (first_sgf, rbs, last_sgf) 02741 CALL array_nullify (row_blk_sizes) 02742 CALL array_new (row_blk_sizes, rbs, gift=.TRUE.) 02743 DEALLOCATE (first_sgf,STAT=istat) 02744 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02745 DEALLOCATE (last_sgf,STAT=istat) 02746 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02747 ! 02748 ! 02749 ALLOCATE(op_mom_ao(1)%matrix) 02750 CALL cp_dbcsr_init(op_mom_ao(1)%matrix,error=error) 02751 CALL cp_dbcsr_create(matrix=op_mom_ao(1)%matrix, & 02752 name="op_mom", & 02753 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry,& 02754 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 02755 nblks=0, nze=0, mutable_work=.TRUE., & 02756 error=error) 02757 CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix,sab_all,error=error) 02758 02759 DO idir2=1,3 02760 ALLOCATE(op_mom_der_ao(1,idir2)%matrix) 02761 CALL cp_dbcsr_init(op_mom_der_ao(1,idir2)%matrix, error=error) 02762 CALL cp_dbcsr_copy(op_mom_der_ao(1,idir2)%matrix,op_mom_ao(1)%matrix,& 02763 "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir2))),error) 02764 ENDDO 02765 02766 DO idir = 2,SIZE(op_mom_ao,1) 02767 ALLOCATE(op_mom_ao(idir)%matrix) 02768 CALL cp_dbcsr_init(op_mom_ao(idir)%matrix, error=error) 02769 CALL cp_dbcsr_copy(op_mom_ao(idir)%matrix,op_mom_ao(1)%matrix,& 02770 "op_mom_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))),error) 02771 DO idir2=1,3 02772 ALLOCATE(op_mom_der_ao(idir,idir2)%matrix) 02773 CALL cp_dbcsr_init(op_mom_der_ao(idir,idir2)%matrix, error=error) 02774 CALL cp_dbcsr_copy(op_mom_der_ao(idir,idir2)%matrix,op_mom_ao(1)%matrix,& 02775 "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir*idir2))),error) 02776 ENDDO 02777 ENDDO 02778 ! 02779 CALL cp_dbcsr_allocate_matrix_set(op_p_ao,3,error=error) 02780 ALLOCATE(op_p_ao(1)%matrix) 02781 CALL cp_dbcsr_init(op_p_ao(1)%matrix,error=error) 02782 CALL cp_dbcsr_create(matrix=op_p_ao(1)%matrix, & 02783 name="op_p_ao", & 02784 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric,& 02785 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 02786 nblks=0, nze=0, mutable_work=.TRUE., & 02787 error=error) 02788 CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb, error=error) 02789 02790 DO idir=2,3 02791 ALLOCATE(op_p_ao(idir)%matrix) 02792 CALL cp_dbcsr_init(op_p_ao(idir)%matrix, error=error) 02793 CALL cp_dbcsr_copy(op_p_ao(idir)%matrix,op_p_ao(1)%matrix,& 02794 "op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))),error) 02795 ENDDO 02796 ! 02797 ! 02798 CALL array_release (row_blk_sizes) 02799 ! 02800 ! 02801 ! Allocate full matrices for only one vector 02802 mo_coeff => psi0_order(1)%matrix 02803 NULLIFY(tmp_fm_struct) 02804 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,& 02805 & ncol_global=max_states,para_env=para_env,& 02806 & context=mo_coeff%matrix_struct%context,& 02807 & error=error) 02808 CALL cp_fm_create(buf,tmp_fm_struct,error=error) 02809 CALL cp_fm_struct_release(tmp_fm_struct,error=error) 02810 ! 02811 ! 02812 ! 02813 ! recompute the linear momentum matrices 02814 CALL build_lin_mom_matrix(qs_env,op_p_ao,error) 02815 !CALL p_xyz_ao(op_p_ao,qs_env,minimum_image=.FALSE.,error=error) 02816 ! 02817 ! 02818 ! get iiB and iiiB 02819 CALL set_vecp(iB,iiB,iiiB) 02820 DO ispin = 1,nspins 02821 ! 02822 CPPostcondition(nbr_center(ispin)==1,cp_failure_level,routineP,error,failure) 02823 ! 02824 ! get ground state MOS 02825 nmo = nstates(ispin) 02826 mo_coeff => psi0_order(ispin)%matrix 02827 CALL get_mo_set(mo_set=mos(ispin)%mo_set,maxocc=maxocc) 02828 ! 02829 ! Initialize the temporary vector chi 02830 chi = 0.0_dp 02831 int_current = 0.0_dp 02832 ! 02833 ! 02834 ! Get the Wannier center of the istate-th ground state orbital 02835 dk(1:3) = centers_set(ispin)%array(1:3,1) 02836 ! 02837 ! Compute the multipole integrals for the state istate, 02838 ! using as reference center the corresponding Wannier center 02839 DO idir = 1,9 02840 CALL cp_dbcsr_set(op_mom_ao(idir)%matrix,0.0_dp,error=error) 02841 DO idir2 = 1,3 02842 CALL cp_dbcsr_set(op_mom_der_ao(idir,idir2)%matrix,0.0_dp,error=error) 02843 ENDDO 02844 ENDDO 02845 CALL rRc_xyz_der_ao(op_mom_ao,op_mom_der_ao,qs_env,dk,order=2,& 02846 & minimum_image=.FALSE.,soft=gapw,error=error) 02847 ! 02848 ! 02849 ! Multuply left and right by the appropriate coefficients and sum into the 02850 ! correct component of the chi tensor using the appropriate multiplicative factor 02851 ! (don't forget the occupation number) 02852 ! Loop over the cartesian components of the tensor 02853 ! The loop over the components of the external field is external, thereby 02854 ! only one column of the chi tensor is computed here 02855 DO idir = 1,3 02856 ! 02857 ! 02858 ! 02859 ! term: dk_ii*2[C0| d_iii(C1(rxp-D))] - dk_iii*2[C0| d_ii(C1(rxp-D))] 02860 IF(.NOT.chi_pbc) THEN 02861 CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix,mo_coeff,& 02862 & buf,ncol=nmo,alpha=1.e0_dp,error=error) 02863 DO jdir = 1,3 02864 DO kdir = 1,3 02865 IF(Levi_Civita(kdir,jdir,idir).EQ.0.0_dp) CYCLE 02866 CALL cp_fm_trace(buf,psi1_rxp(ispin,iB)%matrix,contrib,error=error) 02867 chi(kdir) = chi(kdir) - Levi_Civita(kdir,jdir,idir) * 2.0_dp * dk(jdir) * contrib 02868 ENDDO 02869 ENDDO 02870 ENDIF 02871 ! 02872 ! 02873 ! 02874 ! term: 2[C0| (r-dk)_ii |d_iii(C1(rxp-D))]-2[C0| (r-dk)_iii |d_ii(C1(rxp-D))] 02875 ! and 02876 ! term: -dk_ii * 2[C0|(r-dk)_iiB | d_iii(C1(piiiB))] + 02877 ! +dk_iii * 2[C0|(r-dk)_iiB | d_ii(C1(piiiB))] 02878 ! and 02879 ! term: +dk_ii * 2[C0|(r-dk)_iiiB | d_iii(C1(piiB))] + 02880 ! -dk_iii * 2[C0|(r-dk)_iiiB | d_ii(C1(piiB))] 02881 DO jdir = 1,3 02882 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(jdir,idir)%matrix,mo_coeff,& 02883 & buf,ncol=nmo,alpha=1.e0_dp,error=error) 02884 DO kdir = 1,3 02885 IF(Levi_Civita(kdir,jdir,idir).EQ.0.0_dp) CYCLE 02886 CALL cp_fm_trace(buf,psi1_rxp(ispin,iB)%matrix,contrib,error=error) 02887 chi(kdir) = chi(kdir) - Levi_Civita(kdir,jdir,idir) * 2.0_dp * contrib 02888 ENDDO 02889 ! 02890 IF(.NOT.chi_pbc) THEN 02891 IF(jdir.EQ.iiB) THEN 02892 DO jjdir = 1,3 02893 DO kdir = 1,3 02894 IF(Levi_Civita(kdir,jjdir,idir).EQ.0.0_dp) CYCLE 02895 CALL cp_fm_trace(buf,psi1_p(ispin,iiiB)%matrix,contrib,error=error) 02896 chi(kdir) = chi(kdir) + Levi_Civita(kdir,jjdir,idir) * 2.0_dp * dk(jjdir) * contrib 02897 ENDDO 02898 ENDDO 02899 ENDIF 02900 ! 02901 IF(jdir.EQ.iiiB) THEN 02902 DO jjdir = 1,3 02903 DO kdir = 1,3 02904 IF(Levi_Civita(kdir,jjdir,idir).EQ.0.0_dp) CYCLE 02905 CALL cp_fm_trace(buf,psi1_p(ispin,iiB)%matrix,contrib,error=error) 02906 chi(kdir) = chi(kdir) - Levi_Civita(kdir,jjdir,idir) * 2.0_dp * dk(jjdir) * contrib 02907 ENDDO 02908 ENDDO 02909 ENDIF 02910 ENDIF 02911 ENDDO 02912 ! 02913 ! 02914 ! 02915 ! term1: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] + 02916 ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))] 02917 ! and 02918 ! term1: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] + 02919 ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))] 02920 ! HERE THERE IS ONE EXTRA MULTIPLY 02921 DO jdir = 1,3 02922 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir,iiB),idir)%matrix,mo_coeff,& 02923 & buf,ncol=nmo,alpha=1.e0_dp,error=error) 02924 DO kdir = 1,3 02925 IF(Levi_Civita(kdir,jdir,idir).EQ.0.0_dp) CYCLE 02926 CALL cp_fm_trace(buf,psi1_p(ispin,iiiB)%matrix,contrib,error=error) 02927 chi(kdir) = chi(kdir) + Levi_Civita(kdir,jdir,idir) * 2.0_dp * contrib 02928 ENDDO 02929 ! 02930 CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir,iiiB),idir)%matrix,mo_coeff,& 02931 & buf,ncol=nmo,alpha=1.e0_dp,error=error) 02932 DO kdir = 1,3 02933 IF(Levi_Civita(kdir,jdir,idir).EQ.0.0_dp) CYCLE 02934 CALL cp_fm_trace(buf,psi1_p(ispin,iiB)%matrix,contrib,error=error) 02935 chi(kdir) = chi(kdir) - Levi_Civita(kdir,jdir,idir) * 2.0_dp * contrib 02936 ENDDO 02937 ENDDO 02938 ! 02939 ! 02940 ! 02941 ! term2: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] + 02942 ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))] 02943 ! and 02944 ! term2: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] + 02945 ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))] 02946 CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix,mo_coeff,& 02947 & buf,ncol=nmo,alpha=1.e0_dp,error=error) 02948 DO jdir = 1,3 02949 DO kdir = 1,3 02950 IF(Levi_Civita(kdir,idir,jdir).EQ.0.0_dp) CYCLE 02951 IF(iiB==jdir) THEN 02952 CALL cp_fm_trace(buf,psi1_p(ispin,iiiB)%matrix,contrib,error=error) 02953 chi(kdir) = chi(kdir) + Levi_Civita(kdir,idir,jdir) * contrib 02954 ENDIF 02955 ENDDO 02956 ENDDO 02957 ! 02958 DO jdir = 1,3 02959 DO kdir = 1,3 02960 IF(Levi_Civita(kdir,idir,jdir).EQ.0.0_dp) CYCLE 02961 IF(iiiB==jdir) THEN 02962 CALL cp_fm_trace(buf,psi1_p(ispin,iiB)%matrix,contrib,error=error) 02963 chi(kdir) = chi(kdir) - Levi_Civita(kdir,idir,jdir) * contrib 02964 ENDIF 02965 ! 02966 ENDDO 02967 ENDDO 02968 ! 02969 ! 02970 ! 02971 ! 02972 ENDDO ! idir 02973 ! 02974 DO idir = 1,3 02975 current_env%chi_tensor(idir,iB,ispin) = current_env%chi_tensor(idir,iB,ispin)+& 02976 & maxocc * chi(idir) 02977 IF (output_unit>0) THEN 02978 !WRITE(output_unit,'(A,E12.6)') ' chi_'//ACHAR(119+idir)//ACHAR(119+iB)//& 02979 ! & ' = ',maxocc * chi(idir) 02980 ENDIF 02981 ENDDO 02982 ! 02983 ENDDO ! ispin 02984 ! 02985 ! deallocate the sparse matrices 02986 CALL cp_dbcsr_deallocate_matrix_set(op_mom_ao,error=error) 02987 CALL cp_dbcsr_deallocate_matrix_set(op_mom_der_ao,error=error) 02988 CALL cp_dbcsr_deallocate_matrix_set(op_p_ao,error=error) 02989 CALL cp_fm_release(buf,error=error) 02990 02991 CALL timestop(handle) 02992 02993 END SUBROUTINE current_build_chi_one_center 02994 02995 END MODULE qs_linres_current
1.7.3