|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00011 MODULE qs_linres_epr_ownutils 00012 USE atomic_kind_types, ONLY: atomic_kind_type,& 00013 get_atomic_kind 00014 USE cell_types, ONLY: cell_type 00015 USE cp_control_types, ONLY: dft_control_type 00016 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00017 USE cp_output_handling, ONLY: cp_p_file,& 00018 cp_print_key_finished_output,& 00019 cp_print_key_should_output,& 00020 cp_print_key_unit_nr 00021 USE cp_para_types, ONLY: cp_para_env_type 00022 USE f77_blas 00023 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00024 section_vals_type,& 00025 section_vals_val_get 00026 USE kinds, ONLY: default_string_length,& 00027 dp 00028 USE mathlib, ONLY: diamat_all 00029 USE message_passing, ONLY: mp_sum 00030 USE particle_types, ONLY: particle_type 00031 USE pw_env_types, ONLY: pw_env_get,& 00032 pw_env_type 00033 USE pw_methods, ONLY: pw_axpy,& 00034 pw_integral_ab,& 00035 pw_scale,& 00036 pw_transfer,& 00037 pw_zero 00038 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00039 pw_pool_give_back_pw,& 00040 pw_pool_p_type,& 00041 pw_pool_type 00042 USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc,& 00043 find_coeffs,& 00044 pw_spline_do_precond,& 00045 pw_spline_precond_create,& 00046 pw_spline_precond_release,& 00047 pw_spline_precond_set_kind,& 00048 pw_spline_precond_type,& 00049 spl3_pbc 00050 USE pw_types, ONLY: COMPLEXDATA1D,& 00051 REALDATA3D,& 00052 REALSPACE,& 00053 RECIPROCALSPACE,& 00054 pw_p_type 00055 USE qs_core_energies, ONLY: calculate_ecore 00056 USE qs_environment_types, ONLY: get_qs_env,& 00057 qs_environment_type 00058 USE qs_grid_atom, ONLY: grid_atom_type 00059 USE qs_harmonics_atom, ONLY: harmonics_atom_type 00060 USE qs_linres_nmr_epr_common_utils, ONLY: mult_G_ov_G2_grid 00061 USE qs_linres_op, ONLY: fac_vecp,& 00062 set_vecp,& 00063 set_vecp_rev 00064 USE qs_linres_types, ONLY: current_env_type,& 00065 epr_env_type,& 00066 get_current_env,& 00067 get_epr_env,& 00068 jrho_atom_type,& 00069 nablavks_atom_type 00070 USE qs_rho_atom_types, ONLY: get_rho_atom,& 00071 rho_atom_coeff,& 00072 rho_atom_type 00073 USE qs_rho_types, ONLY: qs_rho_p_type,& 00074 qs_rho_type 00075 USE realspace_grid_types, ONLY: realspace_grid_desc_type 00076 USE timings, ONLY: timeset,& 00077 timestop 00078 USE util, ONLY: get_limit 00079 #include "cp_common_uses.h" 00080 00081 IMPLICIT NONE 00082 00083 PRIVATE 00084 PUBLIC :: epr_g_print, epr_g_zke, epr_g_so, epr_g_soo, epr_ind_magnetic_field 00085 00086 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_ownutils' 00087 00088 CONTAINS 00089 00090 ! ***************************************************************************** 00096 SUBROUTINE epr_g_print(epr_env,qs_env,error) 00097 00098 TYPE(epr_env_type) :: epr_env 00099 TYPE(qs_environment_type), POINTER :: qs_env 00100 TYPE(cp_error_type), INTENT(INOUT), 00101 OPTIONAL :: error 00102 00103 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_print', 00104 routineP = moduleN//':'//routineN 00105 00106 CHARACTER(LEN=default_string_length) :: title 00107 INTEGER :: idir1, idir2, output_unit, 00108 unit_nr 00109 REAL(KIND=dp) :: eigenv_g(3), g_sym(3,3), gsum 00110 TYPE(cp_logger_type), POINTER :: logger 00111 TYPE(section_vals_type), POINTER :: lr_section 00112 00113 NULLIFY(logger, lr_section) 00114 00115 logger => cp_error_get_logger(error) 00116 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00117 00118 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00119 extension=".linresLog",error=error) 00120 00121 gsum = 0.0_dp 00122 00123 DO idir1 = 1,3 00124 DO idir2 = 1,3 00125 gsum = gsum + epr_env%g_total(idir1,idir2) 00126 END DO 00127 END DO 00128 00129 IF (output_unit>0) THEN 00130 WRITE (UNIT=output_unit,FMT="(T2,A,T56,E14.6)")& 00131 "epr|TOT:checksum",gsum 00132 END IF 00133 00134 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 00135 "PRINT%PROGRAM_RUN_INFO",error=error) 00136 00137 IF (BTEST(cp_print_key_should_output(logger%iter_info,lr_section,& 00138 "EPR%PRINT%G_TENSOR",error=error),cp_p_file)) THEN 00139 00140 unit_nr=cp_print_key_unit_nr(logger,lr_section,"EPR%PRINT%G_TENSOR",& 00141 & extension=".data",middle_name="GTENSOR",& 00142 & log_filename=.FALSE.,error=error) 00143 00144 IF(unit_nr > 0) THEN 00145 00146 WRITE(title,"(A)") "G tensor " 00147 WRITE(unit_nr,"(T2,A)") title 00148 00149 WRITE(unit_nr,"(T2,A)") "gmatrix_zke" 00150 WRITE(unit_nr,"(3(A,f15.10))")" XX=",epr_env%g_zke,& 00151 " XY=",0.0_dp," XZ=",0.0_dp 00152 WRITE(unit_nr,"(3(A,f15.10))")" YX=",0.0_dp,& 00153 " YY=",epr_env%g_zke," YZ=",0.0_dp 00154 WRITE(unit_nr,"(3(A,f15.10))")" ZX=",0.0_dp,& 00155 " ZY=",0.0_dp," ZZ=",epr_env%g_zke 00156 00157 WRITE(unit_nr,"(T2,A)") "gmatrix_so" 00158 WRITE(unit_nr,"(3(A,f15.10))")" XX=",epr_env%g_so(1,1),& 00159 " XY=",epr_env%g_so(1,2)," XZ=",epr_env%g_so(1,3) 00160 WRITE(unit_nr,"(3(A,f15.10))")" YX=",epr_env%g_so(2,1),& 00161 " YY=",epr_env%g_so(2,2)," YZ=",epr_env%g_so(2,3) 00162 WRITE(unit_nr,"(3(A,f15.10))")" ZX=",epr_env%g_so(3,1),& 00163 " ZY=",epr_env%g_so(3,2)," ZZ=",epr_env%g_so(3,3) 00164 00165 WRITE(unit_nr,"(T2,A)") "gmatrix_soo" 00166 WRITE(unit_nr,"(3(A,f15.10))")" XX=",epr_env%g_soo(1,1),& 00167 " XY=",epr_env%g_soo(1,2)," XZ=",epr_env%g_soo(1,3) 00168 WRITE(unit_nr,"(3(A,f15.10))")" YX=",epr_env%g_soo(2,1),& 00169 " YY=",epr_env%g_soo(2,2)," YZ=",epr_env%g_soo(2,3) 00170 WRITE(unit_nr,"(3(A,f15.10))")" ZX=",epr_env%g_soo(3,1),& 00171 " ZY=",epr_env%g_soo(3,2)," ZZ=",epr_env%g_soo(3,3) 00172 00173 WRITE(unit_nr,"(T2,A)") "gmatrix_total" 00174 WRITE(unit_nr,"(3(A,f15.10))")" XX=",epr_env%g_total(1,1)+epr_env%g_free_factor,& 00175 " XY=",epr_env%g_total(1,2)," XZ=",epr_env%g_total(1,3) 00176 WRITE(unit_nr,"(3(A,f15.10))")" YX=",epr_env%g_total(2,1),& 00177 " YY=",epr_env%g_total(2,2)+epr_env%g_free_factor," YZ=",epr_env%g_total(2,3) 00178 WRITE(unit_nr,"(3(A,f15.10))")" ZX=",epr_env%g_total(3,1),& 00179 " ZY=",epr_env%g_total(3,2)," ZZ=",epr_env%g_total(3,3)+epr_env%g_free_factor 00180 00181 DO idir1 = 1,3 00182 DO idir2 = 1,3 00183 g_sym(idir1,idir2) = (epr_env%g_total(idir1,idir2) + & 00184 epr_env%g_total(idir2,idir1))/2.0_dp 00185 END DO 00186 END DO 00187 00188 WRITE(unit_nr,"(T2,A)") "gtensor_total" 00189 WRITE(unit_nr,"(3(A,f15.10))")" XX=",g_sym(1,1)+epr_env%g_free_factor,& 00190 " XY=",g_sym(1,2)," XZ=",g_sym(1,3) 00191 WRITE(unit_nr,"(3(A,f15.10))")" YX=",g_sym(2,1),& 00192 " YY=",g_sym(2,2)+epr_env%g_free_factor," YZ=",g_sym(2,3) 00193 WRITE(unit_nr,"(3(A,f15.10))")" ZX=",g_sym(3,1),& 00194 " ZY=",g_sym(3,2)," ZZ=",g_sym(3,3)+epr_env%g_free_factor 00195 00196 CALL diamat_all(g_sym,eigenv_g,error=error) 00197 eigenv_g(:) = eigenv_g(:)*1.0e6_dp 00198 00199 WRITE(unit_nr,"(T2,A)") "delta_g principal values in ppm" 00200 WRITE(unit_nr,"(f15.3,3(A,f15.10))") eigenv_g(1)," X=",g_sym(1,1),& 00201 " Y=",g_sym(2,1)," Z=",g_sym(3,1) 00202 WRITE(unit_nr,"(f15.3,3(A,f15.10))") eigenv_g(2)," X=",g_sym(1,2),& 00203 " Y=",g_sym(2,2)," Z=",g_sym(3,2) 00204 WRITE(unit_nr,"(f15.3,3(A,f15.10))") eigenv_g(3)," X=",g_sym(1,3),& 00205 " Y=",g_sym(2,3)," Z=",g_sym(3,3) 00206 00207 END IF 00208 00209 CALL cp_print_key_finished_output(unit_nr,logger,lr_section,& 00210 & "EPR%PRINT%G_TENSOR",error=error) 00211 00212 ENDIF 00213 00214 END SUBROUTINE epr_g_print 00215 00216 ! ***************************************************************************** 00222 SUBROUTINE epr_g_zke(epr_env,qs_env,error) 00223 00224 TYPE(epr_env_type) :: epr_env 00225 TYPE(qs_environment_type), POINTER :: qs_env 00226 TYPE(cp_error_type), INTENT(INOUT), 00227 OPTIONAL :: error 00228 00229 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_zke', 00230 routineP = moduleN//':'//routineN 00231 00232 INTEGER :: i1, ispin, output_unit 00233 REAL(KIND=dp) :: epr_g_zke_temp(2) 00234 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00235 POINTER :: kinetic 00236 TYPE(cp_logger_type), POINTER :: logger 00237 TYPE(cp_para_env_type), POINTER :: para_env 00238 TYPE(dft_control_type), POINTER :: dft_control 00239 TYPE(qs_rho_type), POINTER :: rho 00240 TYPE(section_vals_type), POINTER :: lr_section 00241 00242 NULLIFY(dft_control, logger, lr_section, rho, kinetic, para_env) 00243 00244 logger => cp_error_get_logger(error) 00245 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00246 00247 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00248 extension=".linresLog",error=error) 00249 00250 CALL get_qs_env(qs_env=qs_env,dft_control=dft_control,& 00251 kinetic=kinetic,rho=rho,para_env=para_env,error=error) 00252 00253 DO ispin=1,dft_control%nspins 00254 CALL calculate_ecore(h=kinetic(1)%matrix,& 00255 p=rho%rho_ao(ispin)%matrix,& 00256 ecore=epr_g_zke_temp(ispin),& 00257 error=error) 00258 END DO 00259 00260 epr_env%g_zke = epr_env%g_zke_factor * ( epr_g_zke_temp(1) - epr_g_zke_temp(2) ) 00261 DO i1 = 1,3 00262 epr_env%g_total(i1,i1) = epr_env%g_total(i1,i1) + epr_env%g_zke 00263 END DO 00264 00265 IF (output_unit>0) THEN 00266 WRITE (UNIT=output_unit,FMT="(T2,A,T56,E24.16)")& 00267 "epr|ZKE:g_zke",epr_env%g_zke 00268 END IF 00269 00270 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 00271 "PRINT%PROGRAM_RUN_INFO",error=error) 00272 00273 END SUBROUTINE epr_g_zke 00274 00275 ! ***************************************************************************** 00281 SUBROUTINE epr_g_so(epr_env,current_env,qs_env,iB,error) 00282 00283 TYPE(epr_env_type) :: epr_env 00284 TYPE(current_env_type) :: current_env 00285 TYPE(qs_environment_type), POINTER :: qs_env 00286 INTEGER, INTENT(IN) :: iB 00287 TYPE(cp_error_type), INTENT(INOUT), 00288 OPTIONAL :: error 00289 00290 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_so', 00291 routineP = moduleN//':'//routineN 00292 00293 INTEGER :: aint_precond, ia, iat, iatom, idir1, idir2, idir3, ikind, ir, 00294 ispin, istat, max_iter, natom, nkind, nspins, output_unit, precond_kind 00295 INTEGER, DIMENSION(2) :: bo 00296 INTEGER, DIMENSION(:), POINTER :: atom_list 00297 LOGICAL :: failure, gapw, paw_atom, 00298 success 00299 REAL(dp) :: eps_r, eps_x, temp_so_soft, 00300 vks_ra_idir2, vks_ra_idir3 00301 REAL(dp), DIMENSION(3, 3) :: temp_so_gapw 00302 REAL(dp), DIMENSION(:, :), POINTER :: g_so, g_total 00303 REAL(KIND=dp), DIMENSION(3) :: ra 00304 TYPE(atomic_kind_type), DIMENSION(:), 00305 POINTER :: atomic_kind_set 00306 TYPE(atomic_kind_type), POINTER :: atom_kind 00307 TYPE(cp_logger_type), POINTER :: logger 00308 TYPE(cp_para_env_type), POINTER :: para_env 00309 TYPE(dft_control_type), POINTER :: dft_control 00310 TYPE(grid_atom_type), POINTER :: grid_atom 00311 TYPE(harmonics_atom_type), POINTER :: harmonics 00312 TYPE(jrho_atom_type), DIMENSION(:), 00313 POINTER :: jrho1_atom_set 00314 TYPE(jrho_atom_type), POINTER :: jrho1_atom 00315 TYPE(nablavks_atom_type), DIMENSION(:), 00316 POINTER :: nablavks_atom_set 00317 TYPE(nablavks_atom_type), POINTER :: nablavks_atom 00318 TYPE(particle_type), DIMENSION(:), 00319 POINTER :: particle_set 00320 TYPE(pw_env_type), POINTER :: pw_env 00321 TYPE(pw_p_type), DIMENSION(:, :), 00322 POINTER :: vks_pw_spline 00323 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00324 TYPE(pw_spline_precond_type), POINTER :: precond 00325 TYPE(qs_rho_p_type), DIMENSION(:), 00326 POINTER :: jrho1_set 00327 TYPE(qs_rho_p_type), DIMENSION(:, :), 00328 POINTER :: nablavks_set 00329 TYPE(section_vals_type), POINTER :: interp_section, lr_section 00330 00331 NULLIFY(atomic_kind_set, atom_kind, atom_list, dft_control,& 00332 grid_atom, g_so, g_total, harmonics, interp_section, jrho1_atom_set,& 00333 jrho1_set, logger, lr_section, nablavks_atom, nablavks_atom_set,& 00334 nablavks_set, para_env, particle_set) 00335 00336 logger => cp_error_get_logger(error) 00337 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00338 00339 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00340 extension=".linresLog",error=error) 00341 00342 failure = .FALSE. 00343 00344 CALL get_qs_env(qs_env=qs_env,dft_control=dft_control,& 00345 atomic_kind_set=atomic_kind_set,& 00346 para_env=para_env,pw_env=pw_env,& 00347 particle_set=particle_set,& 00348 error=error) 00349 00350 CALL get_epr_env(epr_env=epr_env,& 00351 nablavks_set=nablavks_set,& 00352 nablavks_atom_set=nablavks_atom_set,& 00353 g_total=g_total,g_so=g_so) 00354 00355 CALL get_current_env(current_env=current_env,& 00356 jrho1_set=jrho1_set,jrho1_atom_set=jrho1_atom_set,& 00357 error=error) 00358 00359 gapw = dft_control%qs_control%gapw 00360 nkind = SIZE(atomic_kind_set,1) 00361 nspins = dft_control%nspins 00362 00363 DO idir1=1,3 00364 CALL set_vecp(idir1,idir2,idir3) 00365 ! j_pw x nabla_vks_pw 00366 temp_so_soft = 0.0_dp 00367 DO ispin=1,nspins 00368 temp_so_soft = temp_so_soft + (-1.0_dp)**(1 + ispin) * ( & 00369 pw_integral_ab ( jrho1_set(idir2)%rho%rho_r(ispin)%pw, & 00370 nablavks_set(idir3,ispin)%rho%rho_r(1)%pw ,error=error) - & 00371 pw_integral_ab ( jrho1_set(idir3)%rho%rho_r(ispin)%pw, & 00372 nablavks_set(idir2,ispin)%rho%rho_r(1)%pw ,error=error) & 00373 ) 00374 END DO 00375 temp_so_soft = -1.0_dp * epr_env%g_so_factor * temp_so_soft 00376 IF (output_unit>0) THEN 00377 WRITE (UNIT=output_unit,FMT="(T2,A,T18,I1,I1,T56,E24.16)")& 00378 "epr|SOX:soft",iB,idir1,temp_so_soft 00379 END IF 00380 g_so(iB,idir1) = temp_so_soft 00381 END DO !idir1 00382 00383 IF (gapw) THEN 00384 00385 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,error=error) 00386 ALLOCATE(vks_pw_spline(3,nspins),STAT=istat) 00387 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00388 00389 interp_section => section_vals_get_subs_vals(lr_section,& 00390 "EPR%INTERPOLATOR",error=error) 00391 CALL section_vals_val_get(interp_section,"aint_precond", & 00392 i_val=aint_precond, error=error) 00393 CALL section_vals_val_get(interp_section,"precond",i_val=precond_kind,error=error) 00394 CALL section_vals_val_get(interp_section,"max_iter",i_val=max_iter,error=error) 00395 CALL section_vals_val_get(interp_section,"eps_r",r_val=eps_r,error=error) 00396 CALL section_vals_val_get(interp_section,"eps_x",r_val=eps_x,error=error) 00397 00398 DO ispin = 1,nspins 00399 DO idir1 = 1,3 00400 CALL pw_pool_create_pw(auxbas_pw_pool,& 00401 vks_pw_spline(idir1,ispin)%pw,& 00402 use_data=REALDATA3D,in_space=REALSPACE,& 00403 error=error) 00404 ! calculate spline coefficients 00405 CALL pw_spline_precond_create(precond,precond_kind=aint_precond,& 00406 pool=auxbas_pw_pool,pbc=.TRUE.,transpose=.FALSE.,error=error) 00407 CALL pw_spline_do_precond(precond,& 00408 nablavks_set(idir1,ispin)%rho%rho_r(1)%pw,& 00409 vks_pw_spline(idir1,ispin)%pw,error=error) 00410 CALL pw_spline_precond_set_kind(precond,precond_kind,error=error) 00411 success=find_coeffs(values=nablavks_set(idir1,ispin)%rho%rho_r(1)%pw,& 00412 coeffs=vks_pw_spline(idir1,ispin)%pw,linOp=spl3_pbc,& 00413 preconditioner=precond,pool=auxbas_pw_pool,& 00414 eps_r=eps_r,eps_x=eps_x,max_iter=max_iter,& 00415 error=error) 00416 CPPostconditionNoFail(success,cp_warning_level,routineP,error) 00417 CALL pw_spline_precond_release(precond,error=error) 00418 END DO ! idir1 00419 END DO ! ispin 00420 00421 temp_so_gapw = 0.0_dp 00422 00423 DO ikind = 1, nkind 00424 NULLIFY (atom_kind,atom_list,grid_atom,harmonics) 00425 atom_kind => atomic_kind_set(ikind) 00426 CALL get_atomic_kind(atomic_kind=atom_kind,atom_list=atom_list,& 00427 grid_atom=grid_atom,harmonics=harmonics,& 00428 natom=natom, paw_atom=paw_atom) 00429 00430 IF (.NOT.paw_atom) CYCLE 00431 00432 ! Distribute the atoms of this kind 00433 00434 bo = get_limit( natom, para_env%num_pe, para_env%mepos ) 00435 00436 DO iat = 1,natom !bo(1),bo(2)! this partitioning blocks the interpolation 00437 ! ! routines (i.e. waiting for mp_sum) 00438 iatom = atom_list(iat) 00439 NULLIFY (jrho1_atom, nablavks_atom) 00440 jrho1_atom => jrho1_atom_set(iatom) 00441 nablavks_atom => nablavks_atom_set(iatom) 00442 DO idir1 = 1,3 00443 CALL set_vecp(idir1,idir2,idir3) 00444 DO ispin = 1,nspins 00445 DO ir = 1,grid_atom%nr 00446 00447 IF (grid_atom%rad(ir) >= atom_kind%hard_radius) CYCLE 00448 00449 DO ia = 1,grid_atom%ng_sphere 00450 00451 ra = particle_set(iatom)%r 00452 ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:,ia) 00453 vks_ra_idir2 = Eval_Interp_Spl3_pbc(ra,& 00454 vks_pw_spline(idir2,ispin)%pw,error) 00455 vks_ra_idir3 = Eval_Interp_Spl3_pbc(ra,& 00456 vks_pw_spline(idir3,ispin)%pw,error) 00457 00458 IF(iat.LT.bo(1).OR.iat.GT.bo(2))CYCLE!quick and dirty: 00459 ! !here take care of the partition 00460 00461 ! + sum_A j_loc_h_A x nabla_vks_s_A 00462 temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) + & 00463 (-1.0_dp)**(1 + ispin) * ( & 00464 jrho1_atom%jrho_vec_rad_h(idir2,ispin)%r_coef(ir,ia) * & 00465 vks_ra_idir3 - & 00466 jrho1_atom%jrho_vec_rad_h(idir3,ispin)%r_coef(ir,ia) * & 00467 vks_ra_idir2 & 00468 ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00469 00470 ! - sum_A j_loc_s_A x nabla_vks_s_A 00471 temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) - & 00472 (-1.0_dp)**(1 + ispin) * ( & 00473 jrho1_atom%jrho_vec_rad_s(idir2,ispin)%r_coef(ir,ia) * & 00474 vks_ra_idir3 - & 00475 jrho1_atom%jrho_vec_rad_s(idir3,ispin)%r_coef(ir,ia) * & 00476 vks_ra_idir2 & 00477 ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00478 00479 ! + sum_A j_loc_h_A x nabla_vks_loc_h_A 00480 temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) + & 00481 (-1.0_dp)**(1 + ispin) * ( & 00482 jrho1_atom%jrho_vec_rad_h(idir2,ispin)%r_coef(ir,ia) * & 00483 nablavks_atom%nablavks_vec_rad_h(idir3,ispin)%r_coef(ir,ia) - & 00484 jrho1_atom%jrho_vec_rad_h(idir3,ispin)%r_coef(ir,ia) * & 00485 nablavks_atom%nablavks_vec_rad_h(idir2,ispin)%r_coef(ir,ia) & 00486 ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00487 00488 ! - sum_A j_loc_h_A x nabla_vks_loc_s_A 00489 temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) - & 00490 (-1.0_dp)**(1 + ispin) * ( & 00491 jrho1_atom%jrho_vec_rad_h(idir2,ispin)%r_coef(ir,ia) * & 00492 nablavks_atom%nablavks_vec_rad_s(idir3,ispin)%r_coef(ir,ia) - & 00493 jrho1_atom%jrho_vec_rad_h(idir3,ispin)%r_coef(ir,ia) * & 00494 nablavks_atom%nablavks_vec_rad_s(idir2,ispin)%r_coef(ir,ia) & 00495 ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00496 00497 ! ORIGINAL 00498 ! ! + sum_A j_loc_h_A x nabla_vks_loc_h_A 00499 ! temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) + & 00500 ! (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( & 00501 ! jrho1_atom%jrho_vec_rad_h(idir2,iB,ispin)%r_coef(ir,ia) * & 00502 ! nablavks_atom%nablavks_vec_rad_h(idir3,ispin)%r_coef(ir,ia) - & 00503 ! jrho1_atom%jrho_vec_rad_h(idir3,iB,ispin)%r_coef(ir,ia) * & 00504 ! nablavks_atom%nablavks_vec_rad_h(idir2,ispin)%r_coef(ir,ia) & 00505 ! ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00506 ! ! - sum_A j_loc_s_A x nabla_vks_loc_s_A 00507 ! temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) - & 00508 ! (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( & 00509 ! jrho1_atom%jrho_vec_rad_s(idir2,iB,ispin)%r_coef(ir,ia) * & 00510 ! nablavks_atom%nablavks_vec_rad_s(idir3,ispin)%r_coef(ir,ia) - & 00511 ! jrho1_atom%jrho_vec_rad_s(idir3,iB,ispin)%r_coef(ir,ia) * & 00512 ! nablavks_atom%nablavks_vec_rad_s(idir2,ispin)%r_coef(ir,ia) & 00513 ! ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00514 END DO !ia 00515 END DO !ir 00516 END DO !ispin 00517 END DO !idir1 00518 END DO !iat 00519 END DO !ikind 00520 00521 CALL mp_sum(temp_so_gapw,para_env%group) 00522 temp_so_gapw(:,:) = -1.0_dp * epr_env%g_so_factor_gapw * temp_so_gapw(:,:) 00523 00524 IF (output_unit>0) THEN 00525 DO idir1 = 1,3 00526 WRITE (UNIT=output_unit,FMT="(T2,A,T18,I1,I1,T56,E24.16)")& 00527 "epr|SOX:gapw",iB,idir1,temp_so_gapw(iB,idir1) 00528 END DO 00529 END IF 00530 00531 g_so(iB,:) = g_so(iB,:) + temp_so_gapw(iB,:) 00532 00533 DO ispin = 1,nspins 00534 DO idir1 = 1,3 00535 CALL pw_pool_give_back_pw(auxbas_pw_pool,vks_pw_spline(idir1,ispin)%pw,& 00536 error=error) 00537 END DO 00538 END DO 00539 00540 DEALLOCATE(vks_pw_spline,STAT=istat) 00541 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00542 00543 END IF ! gapw 00544 00545 g_total(iB,:) = g_total(iB,:) + g_so(iB,:) 00546 00547 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 00548 "PRINT%PROGRAM_RUN_INFO",error=error) 00549 00550 END SUBROUTINE epr_g_so 00551 00552 ! ***************************************************************************** 00558 SUBROUTINE epr_g_soo(epr_env,current_env,qs_env,iB,error) 00559 00560 TYPE(epr_env_type) :: epr_env 00561 TYPE(current_env_type) :: current_env 00562 TYPE(qs_environment_type), POINTER :: qs_env 00563 INTEGER, INTENT(IN) :: iB 00564 TYPE(cp_error_type), INTENT(INOUT), 00565 OPTIONAL :: error 00566 00567 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_soo', 00568 routineP = moduleN//':'//routineN 00569 00570 INTEGER :: aint_precond, ia, iat, iatom, idir1, ikind, ir, iso, ispin, 00571 istat, max_iter, natom, nkind, nspins, output_unit, precond_kind 00572 INTEGER, DIMENSION(2) :: bo 00573 INTEGER, DIMENSION(:), POINTER :: atom_list 00574 LOGICAL :: failure, gapw, paw_atom, 00575 soo_rho_hard, success 00576 REAL(dp) :: bind_ra_idir1, 00577 chi_tensor(3,3,2), eps_r, 00578 eps_x, rho_spin, temp_soo_soft 00579 REAL(dp), DIMENSION(3, 3) :: temp_soo_gapw 00580 REAL(dp), DIMENSION(:, :), POINTER :: g_soo, g_total 00581 REAL(KIND=dp), DIMENSION(3) :: ra 00582 TYPE(atomic_kind_type), DIMENSION(:), 00583 POINTER :: atomic_kind_set 00584 TYPE(atomic_kind_type), POINTER :: atom_kind 00585 TYPE(cp_logger_type), POINTER :: logger 00586 TYPE(cp_para_env_type), POINTER :: para_env 00587 TYPE(dft_control_type), POINTER :: dft_control 00588 TYPE(grid_atom_type), POINTER :: grid_atom 00589 TYPE(harmonics_atom_type), POINTER :: harmonics 00590 TYPE(particle_type), DIMENSION(:), 00591 POINTER :: particle_set 00592 TYPE(pw_env_type), POINTER :: pw_env 00593 TYPE(pw_p_type), DIMENSION(:, :), 00594 POINTER :: bind_pw_spline 00595 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00596 TYPE(pw_spline_precond_type), POINTER :: precond 00597 TYPE(qs_rho_p_type), DIMENSION(:, :), 00598 POINTER :: bind_set 00599 TYPE(qs_rho_type), POINTER :: rho 00600 TYPE(rho_atom_coeff), DIMENSION(:), 00601 POINTER :: rho_rad_h, rho_rad_s 00602 TYPE(rho_atom_type), DIMENSION(:), 00603 POINTER :: rho_atom_set 00604 TYPE(rho_atom_type), POINTER :: rho_atom 00605 TYPE(section_vals_type), POINTER :: g_section, interp_section, 00606 lr_section 00607 00608 NULLIFY(atomic_kind_set, atom_kind, atom_list, bind_set, dft_control, & 00609 grid_atom, g_section, g_soo, g_total, harmonics, interp_section, & 00610 logger, lr_section, para_env, particle_set, rho, rho_atom, & 00611 rho_atom_set) 00612 00613 logger => cp_error_get_logger(error) 00614 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00615 00616 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00617 extension=".linresLog",error=error) 00618 00619 g_section => section_vals_get_subs_vals(lr_section,& 00620 "EPR%PRINT%G_TENSOR",error=error) 00621 00622 CALL section_vals_val_get(g_section,"soo_rho_hard",l_val=soo_rho_hard,error=error) 00623 00624 failure = .FALSE. 00625 00626 CALL get_qs_env(qs_env=qs_env,atomic_kind_set=atomic_kind_set,& 00627 dft_control=dft_control,para_env=para_env,particle_set=particle_set,& 00628 pw_env=pw_env,rho=rho,rho_atom_set=rho_atom_set,error=error) 00629 00630 CALL get_epr_env(epr_env=epr_env,bind_set=bind_set,& 00631 g_soo=g_soo,g_total=g_total) 00632 00633 CALL get_current_env(current_env=current_env,& 00634 chi_tensor=chi_tensor,& 00635 error=error) 00636 00637 gapw = dft_control%qs_control%gapw 00638 nkind = SIZE(atomic_kind_set,1) 00639 nspins = dft_control%nspins 00640 00641 DO idir1=1,3 00642 temp_soo_soft = 0.0_dp 00643 DO ispin=1,nspins 00644 temp_soo_soft = temp_soo_soft + (-1.0_dp)**(1 +ispin) * ( & 00645 pw_integral_ab ( bind_set(idir1,iB)%rho%rho_r(1)%pw, & 00646 rho%rho_r(ispin)%pw ,error=error) & 00647 ) 00648 END DO 00649 temp_soo_soft = 1.0_dp * epr_env%g_soo_factor * temp_soo_soft 00650 IF (output_unit>0) THEN 00651 WRITE (UNIT=output_unit,FMT="(T2,A,T18,i1,i1,T56,E24.16)")& 00652 "epr|SOO:soft",iB,idir1,temp_soo_soft 00653 END IF 00654 g_soo(iB,idir1) = temp_soo_soft 00655 END DO 00656 00657 DO idir1=1,3 00658 temp_soo_soft = 1.0_dp * epr_env%g_soo_chicorr_factor * chi_tensor(idir1,iB,2) * & 00659 (REAL(dft_control%multiplicity,KIND=dp)-1.0_dp) 00660 IF (output_unit>0) THEN 00661 WRITE (UNIT=output_unit,FMT="(T2,A,T18,i1,i1,T56,E24.16)")& 00662 "epr|SOO:soft_g0",iB,idir1,temp_soo_soft 00663 END IF 00664 g_soo(iB,idir1) = g_soo(iB,idir1) + temp_soo_soft 00665 END DO 00666 00667 IF (gapw .AND. soo_rho_hard) THEN 00668 00669 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,error=error) 00670 ALLOCATE(bind_pw_spline(3,3),STAT=istat) 00671 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00672 00673 interp_section => section_vals_get_subs_vals(lr_section,& 00674 "EPR%INTERPOLATOR",error=error) 00675 CALL section_vals_val_get(interp_section,"aint_precond", & 00676 i_val=aint_precond, error=error) 00677 CALL section_vals_val_get(interp_section,"precond",i_val=precond_kind,error=error) 00678 CALL section_vals_val_get(interp_section,"max_iter",i_val=max_iter,error=error) 00679 CALL section_vals_val_get(interp_section,"eps_r",r_val=eps_r,error=error) 00680 CALL section_vals_val_get(interp_section,"eps_x",r_val=eps_x,error=error) 00681 00682 DO idir1 = 1,3 00683 CALL pw_pool_create_pw(auxbas_pw_pool,& 00684 bind_pw_spline(idir1,iB)%pw,& 00685 use_data=REALDATA3D,in_space=REALSPACE,& 00686 error=error) 00687 ! calculate spline coefficients 00688 CALL pw_spline_precond_create(precond,precond_kind=aint_precond,& 00689 pool=auxbas_pw_pool,pbc=.TRUE.,transpose=.FALSE.,error=error) 00690 CALL pw_spline_do_precond(precond,& 00691 bind_set(idir1,iB)%rho%rho_r(1)%pw,& 00692 bind_pw_spline(idir1,iB)%pw,error=error) 00693 CALL pw_spline_precond_set_kind(precond,precond_kind,error=error) 00694 success=find_coeffs(values=bind_set(idir1,iB)%rho%rho_r(1)%pw,& 00695 coeffs=bind_pw_spline(idir1,iB)%pw,linOp=spl3_pbc,& 00696 preconditioner=precond,pool=auxbas_pw_pool,& 00697 eps_r=eps_r,eps_x=eps_x,max_iter=max_iter,& 00698 error=error) 00699 CPPostconditionNoFail(success,cp_warning_level,routineP,error) 00700 CALL pw_spline_precond_release(precond,error=error) 00701 END DO ! idir1 00702 00703 temp_soo_gapw = 0.0_dp 00704 00705 DO ikind = 1,nkind 00706 NULLIFY (atom_kind,atom_list,grid_atom,harmonics) 00707 atom_kind => atomic_kind_set(ikind) 00708 CALL get_atomic_kind(atomic_kind=atom_kind,atom_list=atom_list,& 00709 grid_atom=grid_atom,harmonics=harmonics,& 00710 natom=natom, paw_atom=paw_atom) 00711 00712 IF (.NOT.paw_atom) CYCLE 00713 00714 ! Distribute the atoms of this kind 00715 00716 bo = get_limit( natom, para_env%num_pe, para_env%mepos ) 00717 00718 DO iat = 1,natom !bo(1),bo(2)! this partitioning blocks the interpolation 00719 ! ! routines (i.e. waiting for mp_sum) 00720 iatom = atom_list(iat) 00721 rho_atom => rho_atom_set(iatom) 00722 NULLIFY(rho_rad_h,rho_rad_s) 00723 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h,& 00724 rho_rad_s=rho_rad_s) 00725 DO idir1 = 1,3 00726 DO ispin = 1,nspins 00727 DO ir = 1,grid_atom%nr 00728 00729 IF (grid_atom%rad(ir) >= atom_kind%hard_radius) CYCLE 00730 00731 DO ia = 1,grid_atom%ng_sphere 00732 00733 ra = particle_set(iatom)%r 00734 ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:,ia) 00735 bind_ra_idir1 = Eval_Interp_Spl3_pbc(ra,& 00736 bind_pw_spline(idir1,iB)%pw,error) 00737 00738 IF(iat.LT.bo(1).OR.iat.GT.bo(2))CYCLE!quick and dirty: 00739 ! !here take care of the partition 00740 00741 rho_spin = 0.0_dp 00742 00743 DO iso = 1,harmonics%max_iso_not0 00744 rho_spin = rho_spin + & 00745 (rho_rad_h(ispin)%r_coef(ir,iso) - & 00746 rho_rad_s(ispin)%r_coef(ir,iso))* & 00747 harmonics%slm(ia,iso) 00748 END DO 00749 00750 temp_soo_gapw(iB,idir1) = temp_soo_gapw(iB,idir1) + & 00751 (-1.0_dp)**(1 + ispin) * ( & 00752 bind_ra_idir1 * rho_spin & 00753 ) * grid_atom%wr(ir)*grid_atom%wa(ia) 00754 00755 END DO !ia 00756 END DO !ir 00757 END DO ! ispin 00758 END DO !idir1 00759 END DO !iat 00760 END DO !ikind 00761 00762 CALL mp_sum(temp_soo_gapw,para_env%group) 00763 temp_soo_gapw(:,:) = 1.0_dp * epr_env%g_soo_factor * temp_soo_gapw(:,:) 00764 00765 IF (output_unit>0) THEN 00766 DO idir1 = 1,3 00767 WRITE (UNIT=output_unit,FMT="(T2,A,T18,I1,I1,T56,E24.16)")& 00768 "epr|SOO:gapw",iB,idir1,temp_soo_gapw(iB,idir1) 00769 END DO 00770 END IF 00771 00772 g_soo(iB,:) = g_soo(iB,:) + temp_soo_gapw(iB,:) 00773 00774 DO idir1 = 1,3 00775 CALL pw_pool_give_back_pw(auxbas_pw_pool,bind_pw_spline(idir1,iB)%pw,& 00776 error=error) 00777 END DO 00778 00779 DEALLOCATE(bind_pw_spline,STAT=istat) 00780 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00781 00782 END IF ! gapw 00783 00784 g_total(iB,:) = g_total(iB,:) + g_soo(iB,:) 00785 00786 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 00787 "PRINT%PROGRAM_RUN_INFO",error=error) 00788 00789 END SUBROUTINE epr_g_soo 00790 00791 SUBROUTINE epr_ind_magnetic_field(epr_env,current_env,qs_env,iB,error) 00792 00793 TYPE(epr_env_type) :: epr_env 00794 TYPE(current_env_type) :: current_env 00795 TYPE(qs_environment_type), POINTER :: qs_env 00796 INTEGER, INTENT(IN) :: iB 00797 TYPE(cp_error_type), INTENT(INOUT) :: error 00798 00799 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_ind_magnetic_field', 00800 routineP = moduleN//':'//routineN 00801 00802 INTEGER :: handle, idir, idir2, idir3, 00803 iiB, iiiB, ispin, istat, 00804 natom, nspins 00805 LOGICAL :: failure, gapw 00806 REAL(dp) :: scale_fac 00807 TYPE(cell_type), POINTER :: cell 00808 TYPE(dft_control_type), POINTER :: dft_control 00809 TYPE(particle_type), DIMENSION(:), 00810 POINTER :: particle_set 00811 TYPE(pw_env_type), POINTER :: pw_env 00812 TYPE(pw_p_type) :: pw_gspace_work, 00813 shift_pw_rspace 00814 TYPE(pw_p_type), DIMENSION(:, :), 00815 POINTER :: shift_pw_gspace 00816 TYPE(pw_p_type), POINTER :: rho_gspace 00817 TYPE(pw_pool_p_type), DIMENSION(:), 00818 POINTER :: pw_pools 00819 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00820 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 00821 00822 failure = .FALSE. 00823 CALL timeset(routineN,handle) 00824 00825 NULLIFY(cell, dft_control, pw_env, auxbas_rs_desc, auxbas_pw_pool, & 00826 rho_gspace, pw_pools, particle_set) 00827 00828 CALL get_qs_env(qs_env=qs_env,cell=cell,dft_control=dft_control,& 00829 particle_set=particle_set,error=error) 00830 00831 gapw = dft_control%qs_control%gapw 00832 natom = SIZE(particle_set,1) 00833 nspins = dft_control%nspins 00834 00835 CALL get_epr_env(epr_env=epr_env,& 00836 error=error) 00837 00838 CALL get_current_env(current_env=current_env,error=error) 00839 00840 00841 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 00842 CALL pw_env_get(pw_env,auxbas_rs_desc=auxbas_rs_desc,& 00843 auxbas_pw_pool=auxbas_pw_pool,pw_pools=pw_pools,& 00844 error=error) 00845 ! 00846 ! Initialize 00847 ! Allocate grids for the calculation of jrho and the shift 00848 ALLOCATE(shift_pw_gspace(3,nspins),STAT=istat) 00849 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00850 DO ispin = 1,nspins 00851 DO idir = 1,3 00852 CALL pw_pool_create_pw(auxbas_pw_pool,shift_pw_gspace(idir,ispin)%pw,& 00853 & use_data=COMPLEXDATA1D,& 00854 & in_space=RECIPROCALSPACE,error=error) 00855 CALL pw_zero(shift_pw_gspace(idir,ispin)%pw,error=error) 00856 ENDDO 00857 ENDDO 00858 CALL pw_pool_create_pw(auxbas_pw_pool,shift_pw_rspace%pw,& 00859 & use_data=REALDATA3D,in_space=REALSPACE,error=error) 00860 CALL pw_zero(shift_pw_rspace%pw,error=error) 00861 CALL pw_pool_create_pw(auxbas_pw_pool,pw_gspace_work%pw,& 00862 & use_data=COMPLEXDATA1D,& 00863 & in_space=RECIPROCALSPACE,error=error) 00864 00865 CALL pw_zero(pw_gspace_work%pw,error=error) 00866 ! 00867 CALL set_vecp(iB,iiB,iiiB) 00868 ! 00869 DO ispin = 1,nspins 00870 ! 00871 DO idir3=1,3 00872 ! set to zero for the calculation of the shift 00873 CALL pw_zero(shift_pw_gspace(idir3,ispin)%pw,error=error) 00874 ENDDO 00875 ! 00876 DO idir = 1,3 00877 rho_gspace => current_env%jrho1_set(idir)%rho%rho_g(ispin) 00878 ! Field gradient 00879 ! loop over the Gvec components: x,y,z 00880 DO idir2 = 1,3 00881 IF(idir /= idir2) THEN 00882 ! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i)) 00883 CALL mult_G_ov_G2_grid(cell,auxbas_pw_pool,rho_gspace,& 00884 & pw_gspace_work,idir2,0.0_dp,error=error) 00885 ! 00886 ! scale and add to the correct component of the shift column 00887 CALL set_vecp_rev(idir,idir2,idir3) 00888 scale_fac=fac_vecp(idir3,idir2,idir) 00889 CALL pw_scale(pw_gspace_work%pw,scale_fac,error=error) 00890 CALL pw_axpy(pw_gspace_work%pw,shift_pw_gspace(idir3,ispin)%pw,error=error) 00891 ENDIF 00892 ENDDO 00893 ! 00894 ENDDO ! idir 00895 ENDDO ! ispin 00896 00897 ! Store the total soft induced magnetic field (corrected for sic) 00898 IF (dft_control%nspins == 2) THEN 00899 DO idir = 1,3 00900 CALL pw_transfer(shift_pw_gspace(idir,2)%pw,epr_env%bind_set(idir,iB)%rho%rho_r(1)%pw,& 00901 error=error) 00902 ENDDO 00903 END IF 00904 ! 00905 ! Dellocate grids for the calculation of jrho and the shift 00906 CALL pw_pool_give_back_pw(auxbas_pw_pool,pw_gspace_work%pw,error=error) 00907 DO ispin = 1,dft_control%nspins 00908 DO idir = 1,3 00909 CALL pw_pool_give_back_pw(auxbas_pw_pool,shift_pw_gspace(idir,ispin)%pw,error=error) 00910 ENDDO 00911 ENDDO 00912 DEALLOCATE(shift_pw_gspace,STAT=istat) 00913 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00914 CALL pw_pool_give_back_pw(auxbas_pw_pool,shift_pw_rspace%pw,error=error) 00915 ! 00916 ! Finalize 00917 CALL timestop(handle) 00918 ! 00919 END SUBROUTINE epr_ind_magnetic_field 00920 00921 END MODULE qs_linres_epr_ownutils 00922
1.7.3