|
CP2K 2.5 (Revision 12981)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00016 MODULE qs_linres_module 00017 USE cp_control_types, ONLY: dft_control_type 00018 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00019 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00020 cp_print_key_unit_nr 00021 USE f77_blas 00022 USE force_env_types, ONLY: force_env_get,& 00023 force_env_type,& 00024 use_qmmm,& 00025 use_qs_force 00026 USE input_constants, ONLY: & 00027 lr_current, lr_none, ot_precond_full_all, ot_precond_full_kinetic, & 00028 ot_precond_full_single, ot_precond_full_single_inverse, & 00029 ot_precond_none, ot_precond_s_inverse, ot_precond_sparse_diag, & 00030 ot_precond_sparse_kinetic, ot_precond_sparse_kinetic_sainv 00031 USE input_section_types, ONLY: section_vals_get,& 00032 section_vals_get_subs_vals,& 00033 section_vals_type,& 00034 section_vals_val_get 00035 USE qmmm_types, ONLY: primary_subsys,& 00036 qs_subsys 00037 USE qs_environment_types, ONLY: get_qs_env,& 00038 qs_environment_type,& 00039 set_qs_env 00040 USE qs_linres_current, ONLY: current_build_chi,& 00041 current_build_current 00042 USE qs_linres_current_utils, ONLY: current_env_cleanup,& 00043 current_env_init,& 00044 current_response 00045 USE qs_linres_epr_nablavks, ONLY: epr_nablavks 00046 USE qs_linres_epr_ownutils, ONLY: epr_g_print,& 00047 epr_g_so,& 00048 epr_g_soo,& 00049 epr_g_zke,& 00050 epr_ind_magnetic_field 00051 USE qs_linres_epr_utils, ONLY: epr_env_cleanup,& 00052 epr_env_init 00053 USE qs_linres_issc_utils, ONLY: issc_env_cleanup,& 00054 issc_env_init,& 00055 issc_issc,& 00056 issc_print,& 00057 issc_response 00058 USE qs_linres_methods, ONLY: linres_localize 00059 USE qs_linres_nmr_shift, ONLY: nmr_shift,& 00060 nmr_shift_print 00061 USE qs_linres_nmr_utils, ONLY: nmr_env_cleanup,& 00062 nmr_env_init 00063 USE qs_linres_op, ONLY: current_operators,& 00064 issc_operators 00065 USE qs_linres_types, ONLY: current_env_type,& 00066 epr_env_type,& 00067 issc_env_type,& 00068 linres_control_create,& 00069 linres_control_release,& 00070 linres_control_type,& 00071 nmr_env_type 00072 USE qs_mo_methods, ONLY: calculate_density_matrix 00073 USE qs_mo_types, ONLY: mo_set_p_type 00074 USE qs_p_env_methods, ONLY: p_env_create,& 00075 p_env_psi0_changed 00076 USE qs_p_env_types, ONLY: p_env_release,& 00077 qs_p_env_type 00078 USE qs_rho_methods, ONLY: qs_rho_update_rho 00079 USE qs_rho_types, ONLY: qs_rho_type 00080 USE termination, ONLY: stop_program 00081 USE timings, ONLY: timeset,& 00082 timestop 00083 #include "cp_common_uses.h" 00084 00085 IMPLICIT NONE 00086 00087 PRIVATE 00088 PUBLIC :: linres_calculation 00089 00090 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_module' 00091 00092 CONTAINS 00093 00094 ! ***************************************************************************** 00103 SUBROUTINE linres_calculation(force_env, error) 00104 00105 TYPE(force_env_type), POINTER :: force_env 00106 TYPE(cp_error_type), INTENT(INOUT) :: error 00107 00108 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation', 00109 routineP = moduleN//':'//routineN 00110 00111 INTEGER :: handle, iatom, iB, output_unit 00112 LOGICAL :: do_qmmm, epr_present, 00113 failure, issc_present, 00114 lr_calculation, nmr_present 00115 TYPE(cp_logger_type), POINTER :: logger 00116 TYPE(current_env_type) :: current_env 00117 TYPE(dft_control_type), POINTER :: dft_control 00118 TYPE(epr_env_type) :: epr_env 00119 TYPE(issc_env_type) :: issc_env 00120 TYPE(linres_control_type), POINTER :: linres_control 00121 TYPE(nmr_env_type) :: nmr_env 00122 TYPE(qs_environment_type), POINTER :: qs_env 00123 TYPE(qs_p_env_type), POINTER :: p_env 00124 TYPE(section_vals_type), POINTER :: lr_section, prop_section 00125 00126 CALL timeset(routineN,handle) 00127 failure = .FALSE. 00128 lr_calculation = .FALSE. 00129 nmr_present = .FALSE. 00130 epr_present = .FALSE. 00131 issc_present = .FALSE. 00132 do_qmmm = .FALSE. 00133 00134 NULLIFY(dft_control,p_env,linres_control,logger,lr_section,prop_section,qs_env) 00135 00136 CPPrecondition(ASSOCIATED(force_env),cp_failure_level,routineP,error,failure) 00137 CPPrecondition(force_env%ref_count>0,cp_failure_level,routineP,error,failure) 00138 00139 SELECT CASE(force_env%in_use) 00140 CASE(use_qs_force) 00141 CALL force_env_get(force_env,& 00142 qs_env=qs_env,error=error) 00143 CASE(use_qmmm) 00144 do_qmmm = .TRUE. 00145 CALL force_env_get(force_env%sub_force_env(primary_subsys)%force_env%sub_force_env(qs_subsys)%force_env,& 00146 qs_env=qs_env,error=error) 00147 CASE DEFAULT 00148 CALL stop_program(routineN,moduleN,__LINE__,"Doesnt recognize this force_env.") 00149 END SELECT 00150 00151 logger => cp_error_get_logger(error) 00152 00153 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00154 CALL section_vals_get(lr_section,explicit=lr_calculation,error=error) 00155 CALL linres_init(lr_section,p_env,qs_env,error=error) 00156 00157 IF(lr_calculation) THEN 00158 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00159 extension=".linresLog",error=error) 00160 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 00161 linres_control=linres_control,error=error) 00162 00163 !The type of perturbation has not been defined yet 00164 linres_control%property = lr_none 00165 ! 00166 ! We do NMR or EPR, then compute the current response 00167 prop_section => section_vals_get_subs_vals(lr_section,"NMR",error=error) 00168 CALL section_vals_get(prop_section,explicit=nmr_present,error=error) 00169 prop_section => section_vals_get_subs_vals(lr_section,"EPR",error=error) 00170 CALL section_vals_get(prop_section,explicit=epr_present,error=error) 00171 00172 IF(nmr_present.OR.epr_present) THEN 00173 linres_control%property = lr_current 00174 IF(.NOT.linres_control%localized_psi0) THEN 00175 CALL stop_program(routineN,moduleN,__LINE__,& 00176 "Are you sure that you want to calculate the chemical "//& 00177 "shift without localized psi0?") 00178 CALL linres_localize(qs_env, linres_control,& 00179 dft_control%nspins,centers_only=.TRUE.,error=error) 00180 ENDIF 00181 IF(dft_control%nspins/=2.AND.epr_present) THEN 00182 CALL stop_program(routineN,moduleN,__LINE__,& 00183 "LSD is needed to perform a g tensor calculation!") 00184 ENDIF 00185 ! 00186 !Initialize the current environment 00187 current_env%ref_count=0 00188 current_env%do_qmmm = do_qmmm 00189 !current_env%prop='nmr' 00190 CALL current_env_init(current_env,qs_env,error=error) 00191 CALL current_operators(current_env,qs_env,error=error) 00192 CALL current_response(current_env,p_env,qs_env,error) 00193 ! 00194 IF(current_env%all_pert_op_done) THEN 00195 !Initialize the nmr environment 00196 IF(nmr_present) THEN 00197 nmr_env%ref_count=0 00198 CALL nmr_env_init(nmr_env,qs_env,error=error) 00199 ENDIF 00200 ! 00201 !Initialize the epr environment 00202 IF(epr_present) THEN 00203 epr_env%ref_count=0 00204 CALL epr_env_init(epr_env,qs_env,error=error) 00205 CALL epr_g_zke(epr_env,qs_env,error=error) 00206 CALL epr_nablavks(epr_env,qs_env,error=error) 00207 ENDIF 00208 ! 00209 ! Build the rs_gauge if needed 00210 !CALL current_set_gauge(current_env,qs_env,error=error) 00211 ! 00212 ! Loop over field direction 00213 DO iB = 1,3 00214 ! 00215 ! Build current response and succeptibility 00216 CALL current_build_current(current_env,qs_env,iB,error=error) 00217 CALL current_build_chi(current_env,qs_env,iB,error=error) 00218 ! 00219 ! Compute NMR shift 00220 IF(nmr_present) THEN 00221 CALL nmr_shift(nmr_env,current_env,qs_env,iB,error=error) 00222 ENDIF 00223 ! 00224 ! Compute EPR 00225 IF(epr_present) THEN 00226 CALL epr_ind_magnetic_field(epr_env,current_env,qs_env,iB,error=error) 00227 CALL epr_g_so(epr_env,current_env,qs_env,iB,error=error) 00228 CALL epr_g_soo(epr_env,current_env,qs_env,iB,error=error) 00229 ENDIF 00230 ENDDO 00231 ! 00232 ! Finalized the nmr environment 00233 IF(nmr_present) THEN 00234 CALL nmr_shift_print(nmr_env,current_env,qs_env,error=error) 00235 CALL nmr_env_cleanup(nmr_env,error=error) 00236 ENDIF 00237 ! 00238 ! Finalized the epr environment 00239 IF(epr_present) THEN 00240 CALL epr_g_print(epr_env,qs_env,error=error) 00241 CALL epr_env_cleanup(epr_env,error=error) 00242 ENDIF 00243 ! 00244 ELSE 00245 IF(output_unit>0) THEN 00246 WRITE(output_unit,"(T10,A,/T20,A,/)")& 00247 "CURRENT: Not all responses to perturbation operators could be calculated.",& 00248 " Hence: NO nmr and NO epr possible." 00249 END IF 00250 END IF 00251 ! Finalized the current environment 00252 CALL current_env_cleanup(current_env,qs_env,error=error) 00253 ENDIF 00254 ! 00255 ! We do the indirect spin-spin coupling calculation 00256 prop_section => section_vals_get_subs_vals(lr_section,"SPINSPIN",error=error) 00257 CALL section_vals_get(prop_section,explicit=issc_present,error=error) 00258 00259 IF(issc_present) THEN 00260 linres_control%property = lr_current 00261 IF(.NOT.linres_control%localized_psi0) THEN 00262 CALL stop_program(routineN,moduleN,__LINE__,& 00263 "Are you sure that you want to calculate the chemical "//& 00264 "shift without localized psi0?") 00265 CALL linres_localize(qs_env,linres_control,& 00266 dft_control%nspins,centers_only=.TRUE.,error=error) 00267 ENDIF 00268 ! 00269 !Initialize the current environment 00270 current_env%ref_count=0 00271 current_env%do_qmmm = do_qmmm 00272 !current_env%prop='issc' 00273 !CALL current_env_init(current_env,qs_env,error=error) 00274 !CALL current_response(current_env,p_env,qs_env,error) 00275 ! 00276 !Initialize the issc environment 00277 issc_env%ref_count=0 00278 CALL issc_env_init(issc_env,qs_env,error=error) 00279 ! 00280 ! Loop over atoms 00281 DO iatom = 1,issc_env%issc_natms 00282 CALL issc_operators(issc_env,qs_env,iatom,error) 00283 CALL issc_response(issc_env,p_env,qs_env,error) 00284 CALL issc_issc(issc_env,qs_env,iatom,error=error) 00285 ENDDO 00286 ! 00287 ! Finalized the issc environment 00288 CALL issc_print(issc_env,qs_env,error) 00289 CALL issc_env_cleanup(issc_env,error) 00290 ENDIF 00291 00292 ! Other possible LR calculations can be introduced here 00293 00294 CALL p_env_release(p_env,error=error) 00295 00296 IF(output_unit>0) THEN 00297 WRITE (UNIT=output_unit,FMT="(/,T3,A,/,T25,A,/,T3,A,/)")& 00298 REPEAT("=",77),& 00299 "ENDED LINRES CALCULATION",& 00300 REPEAT("=",77) 00301 END IF 00302 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 00303 "PRINT%PROGRAM_RUN_INFO",error=error) 00304 ELSE 00305 output_unit = cp_logger_get_default_io_unit(logger) 00306 IF(output_unit>0) THEN 00307 WRITE (output_unit, "(2X,A)") "",& 00308 "-----------------------------------------------------------------------------",& 00309 "- No LR calculation has been specified in the input -",& 00310 " cp2k is going to stop, bye bye ",& 00311 "-----------------------------------------------------------------------------",& 00312 "" 00313 END IF 00314 END IF 00315 00316 CALL timestop(handle) 00317 00318 END SUBROUTINE linres_calculation 00319 00320 ! ***************************************************************************** 00329 SUBROUTINE linres_init(lr_section,p_env,qs_env,error) 00330 00331 TYPE(section_vals_type), POINTER :: lr_section 00332 TYPE(qs_p_env_type), POINTER :: p_env 00333 TYPE(qs_environment_type), POINTER :: qs_env 00334 TYPE(cp_error_type), INTENT(INOUT) :: error 00335 00336 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_init', 00337 routineP = moduleN//':'//routineN 00338 00339 INTEGER :: ispin, output_unit 00340 LOGICAL :: do_it, failure 00341 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00342 POINTER :: matrix_ks 00343 TYPE(cp_logger_type), POINTER :: logger 00344 TYPE(dft_control_type), POINTER :: dft_control 00345 TYPE(linres_control_type), POINTER :: linres_control 00346 TYPE(mo_set_p_type), DIMENSION(:), 00347 POINTER :: mos 00348 TYPE(qs_rho_type), POINTER :: rho 00349 TYPE(section_vals_type), POINTER :: loc_section 00350 00351 failure = .FALSE. 00352 00353 NULLIFY(logger) 00354 logger => cp_error_get_logger(error) 00355 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00356 extension=".linresLog",error=error) 00357 NULLIFY(dft_control, linres_control, loc_section, rho, mos, matrix_ks) 00358 00359 CPPrecondition(.NOT.ASSOCIATED(p_env),cp_failure_level,routineP,error,failure) 00360 IF (.NOT. failure) THEN 00361 00362 CALL linres_control_create(linres_control,error=error) 00363 CALL set_qs_env(qs_env=qs_env, linres_control=linres_control,error=error) 00364 CALL linres_control_release(linres_control,error=error) 00365 CALL get_qs_env(qs_env=qs_env, linres_control=linres_control,& 00366 dft_control=dft_control,matrix_ks=matrix_ks,mos=mos,rho=rho,error=error) 00367 00368 ! Localized Psi0 are required when the position operator has to be defined (nmr) 00369 loc_section =>section_vals_get_subs_vals(lr_section,"LOCALIZE",error=error) 00370 CALL section_vals_val_get(loc_section,"_SECTION_PARAMETERS_",& 00371 l_val=linres_control%localized_psi0,error=error) 00372 IF(linres_control%localized_psi0) THEN 00373 IF(output_unit>0) THEN 00374 WRITE (UNIT=output_unit,FMT="(/,T3,A,A)")& 00375 "Localization of the ground state orbitals",& 00376 " before starting the linear response calculation" 00377 END IF 00378 00379 CALL linres_localize(qs_env, linres_control,dft_control%nspins,error=error) 00380 00381 DO ispin=1,dft_control%nspins 00382 CALL calculate_density_matrix(mos(ispin)%mo_set,rho%rho_ao(ispin)%matrix,error=error) 00383 ENDDO 00384 ! ** update qs_env%rho 00385 CALL qs_rho_update_rho(rho, qs_env=qs_env, error=error) 00386 END IF 00387 00388 CALL section_vals_val_get(lr_section,"RESTART",l_val=linres_control%linres_restart,error=error) 00389 CALL section_vals_val_get(lr_section,"MAX_ITER",i_val=linres_control%max_iter,error=error) 00390 CALL section_vals_val_get(lr_section,"EPS",r_val=linres_control%eps,error=error) 00391 CALL section_vals_val_get(lr_section,"RESTART_EVERY",i_val=linres_control%restart_every,error=error) 00392 CALL section_vals_val_get(lr_section,"PRECONDITIONER",i_val=linres_control%preconditioner_type,error=error) 00393 CALL section_vals_val_get(lr_section,"ENERGY_GAP",r_val=linres_control%energy_gap,error=error) 00394 00395 IF(output_unit>0) THEN 00396 WRITE (UNIT=output_unit,FMT="(/,T3,A,/,T25,A,/,T3,A,/)")& 00397 REPEAT("=",77),& 00398 "START LINRES CALCULATION",& 00399 REPEAT("=",77) 00400 00401 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 00402 "Properties to be Calulated:" 00403 CALL section_vals_val_get(lr_section,"NMR%_SECTION_PARAMETERS_",& 00404 l_val=do_it,error=error) 00405 IF(do_it) WRITE (UNIT=output_unit,FMT="(T45,A)") & 00406 "NMR Chemical Shift" 00407 00408 IF(linres_control%localized_psi0) WRITE (UNIT=output_unit,FMT="(T2,A,T65,A)")& 00409 "LINRES|"," LOCALIZED PSI0" 00410 00411 WRITE(UNIT=output_unit,FMT="(T2,A,T60,A)")& 00412 "LINRES| Optimization algorithm"," Conjugate Gradients" 00413 00414 SELECT CASE(linres_control%preconditioner_type) 00415 CASE(ot_precond_none) 00416 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00417 "LINRES| Preconditioner"," NONE" 00418 CASE(ot_precond_full_single) 00419 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00420 "LINRES| Preconditioner"," FULL_SINGLE" 00421 CASE(ot_precond_full_kinetic) 00422 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00423 "LINRES| Preconditioner"," FULL_KINETIC" 00424 CASE(ot_precond_s_inverse) 00425 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00426 "LINRES| Preconditioner"," FULL_S_INVERSE" 00427 CASE(ot_precond_full_single_inverse) 00428 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00429 "LINRES| Preconditioner"," FULL_SINGLE_INVERSE" 00430 CASE(ot_precond_full_all) 00431 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00432 "LINRES| Preconditioner"," FULL_ALL" 00433 CASE(ot_precond_sparse_diag) 00434 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00435 "LINRES| Preconditioner"," SPARSE_DIAG" 00436 CASE(ot_precond_sparse_kinetic) 00437 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00438 "LINRES| Preconditioner"," SPARSE_KINETIC" 00439 CASE(ot_precond_sparse_kinetic_sainv) 00440 WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")& 00441 "LINRES| Preconditioner","SAINV SPARSE KINETIC" 00442 CASE DEFAULT 00443 CALL stop_program(routineN,moduleN,__LINE__,"Preconditioner NYI") 00444 END SELECT 00445 00446 WRITE (UNIT=output_unit,FMT="(T2,A,T72,ES8.1)")& 00447 "LINRES| EPS",linres_control%eps 00448 WRITE (UNIT=output_unit,FMT="(T2,A,T72,I8)")& 00449 "LINRES| MAX_ITER",linres_control%max_iter 00450 END IF 00451 00452 !------------------! 00453 ! create the p_env ! 00454 !------------------! 00455 CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.TRUE.,linres_control=linres_control,error=error) 00456 00457 ! update the m_epsilon matrix 00458 CALL p_env_psi0_changed(p_env,qs_env,error=error) 00459 00460 ! calculate eigenvectros and eigenvalues of K 00461 p_env%os_valid = .FALSE. 00462 p_env%new_preconditioner = .TRUE. 00463 END IF 00464 CALL cp_print_key_finished_output(output_unit,logger,lr_section,& 00465 "PRINT%PROGRAM_RUN_INFO",error=error) 00466 00467 END SUBROUTINE linres_init 00468 00469 END MODULE qs_linres_module 00470
1.7.3