|
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 ! ***************************************************************************** 00013 MODULE qs_linres_methods 00014 USE cell_types, ONLY: cell_type 00015 USE cp_control_types, ONLY: dft_control_type 00016 USE cp_dbcsr_interface, ONLY: cp_dbcsr_checksum,& 00017 cp_dbcsr_copy,& 00018 cp_dbcsr_init,& 00019 cp_dbcsr_set 00020 USE cp_dbcsr_operations, ONLY: cp_dbcsr_allocate_matrix_set,& 00021 cp_dbcsr_plus_fm_fm_t,& 00022 cp_dbcsr_sm_fm_multiply 00023 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00024 cp_dbcsr_type 00025 USE cp_external_control, ONLY: external_control 00026 USE cp_files, ONLY: close_file,& 00027 open_file 00028 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm,& 00029 cp_fm_scale_and_add,& 00030 cp_fm_trace 00031 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00032 cp_fm_struct_release,& 00033 cp_fm_struct_type 00034 USE cp_fm_types, ONLY: cp_fm_create,& 00035 cp_fm_get_info,& 00036 cp_fm_get_submatrix,& 00037 cp_fm_p_type,& 00038 cp_fm_release,& 00039 cp_fm_set_submatrix,& 00040 cp_fm_to_fm,& 00041 cp_fm_type 00042 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 00043 USE cp_output_handling, ONLY: cp_p_file,& 00044 cp_print_key_finished_output,& 00045 cp_print_key_generate_filename,& 00046 cp_print_key_should_output,& 00047 cp_print_key_unit_nr 00048 USE cp_para_types, ONLY: cp_para_env_type 00049 USE f77_blas 00050 USE hartree_local_methods, ONLY: Vh_1c_gg_integrals 00051 USE input_constants, ONLY: do_loc_none,& 00052 op_loc_berry,& 00053 ot_precond_none,& 00054 ot_precond_solver_default,& 00055 state_loc_all 00056 USE input_section_types, ONLY: section_get_ival,& 00057 section_get_rval,& 00058 section_vals_get_subs_vals,& 00059 section_vals_type,& 00060 section_vals_val_get 00061 USE kinds, ONLY: default_path_length,& 00062 default_string_length,& 00063 dp 00064 USE machine, ONLY: m_flush,& 00065 m_walltime 00066 USE message_passing, ONLY: mp_bcast 00067 USE preconditioner, ONLY: apply_preconditioner,& 00068 make_preconditioner 00069 USE pw_env_types, ONLY: pw_env_get,& 00070 pw_env_type 00071 USE pw_methods, ONLY: pw_axpy,& 00072 pw_copy,& 00073 pw_transfer,& 00074 pw_zero 00075 USE pw_poisson_methods, ONLY: pw_poisson_solve 00076 USE pw_poisson_types, ONLY: pw_poisson_type 00077 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00078 pw_pool_give_back_pw,& 00079 pw_pool_p_type,& 00080 pw_pool_type 00081 USE pw_types, ONLY: COMPLEXDATA1D,& 00082 REALDATA3D,& 00083 REALSPACE,& 00084 RECIPROCALSPACE,& 00085 pw_create,& 00086 pw_p_type,& 00087 pw_release,& 00088 pw_retain 00089 USE qs_environment_types, ONLY: get_qs_env,& 00090 qs_environment_type 00091 USE qs_gapw_densities, ONLY: prepare_gapw_den 00092 USE qs_integrate_potential, ONLY: integrate_v_rspace 00093 USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type 00094 USE qs_ks_atom, ONLY: update_ks_atom 00095 USE qs_linres_types, ONLY: linres_control_type 00096 USE qs_loc_methods, ONLY: qs_loc_driver 00097 USE qs_loc_types, ONLY: get_qs_loc_env,& 00098 localized_wfn_control_type,& 00099 qs_loc_env_create,& 00100 qs_loc_env_new_type,& 00101 qs_loc_env_release,& 00102 qs_loc_env_retain 00103 USE qs_loc_utils, ONLY: loc_write_restart,& 00104 qs_loc_control_init,& 00105 qs_loc_init 00106 USE qs_mo_types, ONLY: get_mo_set,& 00107 mo_set_p_type 00108 USE qs_p_env_types, ONLY: qs_p_env_type 00109 USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace 00110 USE qs_rho_methods, ONLY: qs_rho_rebuild,& 00111 qs_rho_update_rho 00112 USE qs_rho_types, ONLY: qs_rho_get,& 00113 qs_rho_type 00114 USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom 00115 USE string_utilities, ONLY: xstring 00116 USE termination, ONLY: stop_program 00117 USE timings, ONLY: timeset,& 00118 timestop 00119 USE xc, ONLY: xc_calc_2nd_deriv,& 00120 xc_prep_2nd_deriv 00121 USE xc_derivatives, ONLY: xc_functionals_get_needs 00122 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00123 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 00124 xc_rho_set_release,& 00125 xc_rho_set_type,& 00126 xc_rho_set_update 00127 #include "cp_common_uses.h" 00128 00129 IMPLICIT NONE 00130 00131 PRIVATE 00132 00133 ! *** Public subroutines *** 00134 PUBLIC :: linres_localize, linres_solver 00135 PUBLIC :: linres_write_restart, linres_read_restart, linres_init_write_restart 00136 00137 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_methods' 00138 00139 CONTAINS 00140 00141 ! ***************************************************************************** 00148 SUBROUTINE linres_localize(qs_env,linres_control,nspins,centers_only,error) 00149 00150 TYPE(qs_environment_type), POINTER :: qs_env 00151 TYPE(linres_control_type), POINTER :: linres_control 00152 INTEGER, INTENT(IN) :: nspins 00153 LOGICAL, INTENT(IN), OPTIONAL :: centers_only 00154 TYPE(cp_error_type), INTENT(INOUT) :: error 00155 00156 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_localize', 00157 routineP = moduleN//':'//routineN 00158 00159 INTEGER :: ispin, istat, istate, 00160 nmoloc(2), output_unit 00161 LOGICAL :: failure, my_centers_only 00162 TYPE(cp_fm_p_type), DIMENSION(:), 00163 POINTER :: mos_localized 00164 TYPE(cp_fm_type), POINTER :: mo_coeff 00165 TYPE(cp_logger_type), POINTER :: logger 00166 TYPE(localized_wfn_control_type), 00167 POINTER :: localized_wfn_control 00168 TYPE(mo_set_p_type), DIMENSION(:), 00169 POINTER :: mos 00170 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 00171 TYPE(section_vals_type), POINTER :: loc_print_section, 00172 loc_section, lr_section 00173 00174 failure = .FALSE. 00175 NULLIFY(logger, lr_section, loc_section, loc_print_section,localized_wfn_control) 00176 logger => cp_error_get_logger(error) 00177 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00178 loc_section => section_vals_get_subs_vals(lr_section,"LOCALIZE",error=error) 00179 loc_print_section => section_vals_get_subs_vals(lr_section,"LOCALIZE%PRINT",error=error) 00180 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00181 extension=".linresLog",error=error) 00182 my_centers_only = .FALSE. 00183 IF(PRESENT(centers_only)) my_centers_only = centers_only 00184 00185 NULLIFY(mos,mo_coeff,qs_loc_env,mos_localized) 00186 CALL get_qs_env(qs_env=qs_env,mos=mos,error=error) 00187 CALL qs_loc_env_create(qs_loc_env,error=error) 00188 CALL qs_loc_env_retain(qs_loc_env, error=error) 00189 linres_control% qs_loc_env=> qs_loc_env 00190 CALL qs_loc_env_release(qs_loc_env,error=error) 00191 qs_loc_env => linres_control% qs_loc_env 00192 CALL qs_loc_control_init(qs_loc_env,loc_section,do_homo=.TRUE.,error=error) 00193 CALL get_qs_loc_env(qs_loc_env,localized_wfn_control=localized_wfn_control,error=error) 00194 00195 ALLOCATE(mos_localized(nspins),stat=istat) 00196 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00197 DO ispin = 1,nspins 00198 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff) 00199 CALL cp_fm_create(mos_localized(ispin)%matrix,mo_coeff%matrix_struct,error=error) 00200 CALL cp_fm_to_fm(mo_coeff,mos_localized(ispin)%matrix,error=error) 00201 END DO 00202 00203 00204 nmoloc(1:2) = 0 00205 IF(my_centers_only) THEN 00206 localized_wfn_control%set_of_states = state_loc_all 00207 localized_wfn_control%localization_method = do_loc_none 00208 localized_wfn_control%operator_type = op_loc_berry 00209 ENDIF 00210 00211 CALL qs_loc_init(qs_env, qs_loc_env,loc_section,mos_localized=mos_localized,& 00212 do_homo=.TRUE.,error=error) 00213 00214 ! The orbital centers are stored in linres_control%localized_wfn_control 00215 DO ispin = 1,nspins 00216 CALL qs_loc_driver(qs_env,qs_loc_env,loc_section,loc_print_section, myspin = ispin,& 00217 ext_mo_coeff=mos_localized(ispin)%matrix,error=error) 00218 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff) 00219 CALL cp_fm_to_fm(mos_localized(ispin)%matrix,mo_coeff,error=error) 00220 END DO 00221 00222 CALL loc_write_restart(qs_env,qs_loc_env,loc_print_section,mos,& 00223 mos_localized, do_homo=.TRUE., error=error) 00224 CALL cp_fm_vect_dealloc(mos_localized,error) 00225 00226 00227 ! Write Centers and Spreads on std out 00228 IF(output_unit > 0) THEN 00229 DO ispin = 1,nspins 00230 WRITE (output_unit,"(/,T2,A,I2)")& 00231 "WANNIER CENTERS for spin ",ispin 00232 WRITE (output_unit,"(/,T18,A,3X,A)")& 00233 "--------------- Centers --------------- ",& 00234 "--- Spreads ---" 00235 DO istate = 1,SIZE(localized_wfn_control%centers_set(ispin)%array,2) 00236 WRITE(output_unit,"(T5,A6,I6,2X,3f12.6,5X,f12.6)")& 00237 'state ', istate,localized_wfn_control%centers_set(ispin)%array(1:3,istate),& 00238 localized_wfn_control%centers_set(ispin)%array(4,istate) 00239 END DO 00240 END DO 00241 CALL m_flush(output_unit) 00242 END IF 00243 00244 00245 END SUBROUTINE linres_localize 00246 00247 ! ***************************************************************************** 00255 SUBROUTINE linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, should_stop, error) 00256 ! 00257 TYPE(qs_p_env_type), POINTER :: p_env 00258 TYPE(qs_environment_type), POINTER :: qs_env 00259 TYPE(cp_fm_p_type), DIMENSION(:), 00260 POINTER :: psi1, h1_psi0, psi0_order 00261 LOGICAL, INTENT(OUT) :: should_stop 00262 TYPE(cp_error_type), INTENT(INOUT) :: error 00263 00264 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_solver', 00265 routineP = moduleN//':'//routineN 00266 00267 INTEGER :: handle, ispin, istat, iter, 00268 maxnmo, maxnmo_o, nao, ncol, 00269 nmo, nspins, output_unit 00270 LOGICAL :: failure, restart 00271 REAL(dp) :: norm_res, t1, t2 00272 REAL(dp), DIMENSION(:), POINTER :: alpha, beta, tr_pAp, tr_rz0, 00273 tr_rz1 00274 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00275 POINTER :: matrix_ks, matrix_s, matrix_t 00276 TYPE(cp_fm_p_type), DIMENSION(:), 00277 POINTER :: Ap, chc, mo_coeff_array, p, 00278 r, Sc, z 00279 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 00280 TYPE(cp_fm_type), POINTER :: buf, mo_coeff 00281 TYPE(cp_logger_type), POINTER :: logger 00282 TYPE(cp_para_env_type), POINTER :: para_env 00283 TYPE(dft_control_type), POINTER :: dft_control 00284 TYPE(linres_control_type), POINTER :: linres_control 00285 TYPE(mo_set_p_type), DIMENSION(:), 00286 POINTER :: mos 00287 TYPE(section_vals_type), POINTER :: lr_section 00288 00289 ! 00290 00291 failure = .FALSE. 00292 ! 00293 CALL timeset(routineN,handle) 00294 00295 NULLIFY(dft_control,linres_control,matrix_s,matrix_t,matrix_ks,para_env) 00296 NULLIFY(Ap,r,p,z,lr_section,logger,buf,mos,tmp_fm_struct,mo_coeff) 00297 NULLIFY(Sc,chc) 00298 00299 logger => cp_error_get_logger(error) 00300 t1 = m_walltime() 00301 00302 CALL get_qs_env(qs_env=qs_env,& 00303 matrix_ks=matrix_ks,& 00304 matrix_s=matrix_s,& 00305 kinetic=matrix_t,& 00306 dft_control=dft_control,& 00307 linres_control=linres_control,& 00308 para_env=para_env,& 00309 mos=mos,& 00310 error=error) 00311 00312 ! 00313 nspins = dft_control%nspins 00314 CALL get_mo_set(mos(1)%mo_set,nao=nao) 00315 maxnmo = 0 00316 maxnmo_o = 0 00317 DO ispin = 1,nspins 00318 CALL get_mo_set(mos(ispin)%mo_set,nmo=ncol) 00319 maxnmo = MAX(maxnmo,ncol) 00320 CALL cp_fm_get_info(psi0_order(ispin)%matrix,ncol_global=ncol,error=error) 00321 maxnmo_o = MAX(maxnmo_o,ncol) 00322 ENDDO 00323 ! 00324 lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error) 00325 output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",& 00326 extension=".linresLog",error=error) 00327 00328 ! 00329 CALL check_p_env_init(p_env,linres_control,nspins,error=error) 00330 ! 00331 ! allocate the vectors 00332 ALLOCATE(alpha(nspins),beta(nspins),tr_pAp(nspins),tr_rz0(nspins),tr_rz1(nspins),& 00333 r(nspins),p(nspins),z(nspins),Ap(nspins),mo_coeff_array(nspins),STAT=istat) 00334 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00335 DO ispin = 1,nspins 00336 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 00337 mo_coeff_array(ispin)%matrix => mo_coeff 00338 ENDDO 00339 ! 00340 DO ispin = 1,nspins 00341 NULLIFY(r(ispin)%matrix,p(ispin)%matrix,z(ispin)%matrix,Ap(ispin)%matrix) 00342 CALL cp_fm_create(r(ispin)%matrix,psi1(ispin)%matrix%matrix_struct,error=error) 00343 CALL cp_fm_create(p(ispin)%matrix,psi1(ispin)%matrix%matrix_struct,error=error) 00344 CALL cp_fm_create(z(ispin)%matrix,psi1(ispin)%matrix%matrix_struct,error=error) 00345 CALL cp_fm_create(Ap(ispin)%matrix,psi1(ispin)%matrix%matrix_struct,error=error) 00346 ENDDO 00347 ! 00348 NULLIFY(tmp_fm_struct) 00349 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,& 00350 & ncol_global=maxnmo,para_env=para_env,& 00351 & context=psi1(1)%matrix%matrix_struct%context,& 00352 & error=error) 00353 CALL cp_fm_create(buf,tmp_fm_struct,error=error) 00354 CALL cp_fm_struct_release(tmp_fm_struct,error=error) 00355 ! 00356 ! 00357 ! 00358 ! compute S*C0, C0_order'*H*C0_order (this should be done once for all) 00359 ALLOCATE(chc(nspins),Sc(nspins),STAT=istat) 00360 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00361 DO ispin = 1,nspins 00362 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff,nmo=nmo) 00363 CALL cp_fm_create(Sc(ispin)%matrix,mo_coeff%matrix_struct,error=error) 00364 NULLIFY(tmp_fm_struct) 00365 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nmo,& 00366 & ncol_global=nmo,para_env=para_env,& 00367 & context=mo_coeff%matrix_struct%context,& 00368 & error=error) 00369 CALL cp_fm_create(chc(ispin)%matrix,tmp_fm_struct,error=error) 00370 CALL cp_fm_struct_release(tmp_fm_struct,error=error) 00371 ENDDO 00372 ! 00373 DO ispin = 1,nspins 00374 ! 00375 ! C0_order' * H * C0_order 00376 mo_coeff => psi0_order(ispin)%matrix 00377 CALL cp_fm_get_info(mo_coeff,ncol_global=ncol,error=error) 00378 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix,mo_coeff,buf,ncol,error=error) 00379 CALL cp_fm_gemm('T','N',ncol,ncol,nao,-1.0_dp,mo_coeff,buf,0.0_dp,chc(ispin)%matrix,error) 00380 ! 00381 ! S * C0 00382 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 00383 CALL cp_fm_get_info(mo_coeff,ncol_global=ncol,error=error) 00384 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix,mo_coeff,Sc(ispin)%matrix,ncol,error=error) 00385 ENDDO 00386 ! 00387 ! 00388 ! 00389 ! header 00390 IF(output_unit>0) THEN 00391 WRITE(output_unit,"(/,T3,A,T16,A,T25,A,T38,A,T52,A,T72,A,/,T3,A)")& 00392 "Iteration","Method","Restart","Stepsize","Convergence","Time",& 00393 REPEAT("-",80) 00394 ENDIF 00395 ! 00396 ! orthogonalize x with respect to the psi0 00397 CALL preortho(psi1,mo_coeff_array,Sc,buf,error) 00398 ! 00399 ! build the preconditioner 00400 IF(linres_control%preconditioner_type /= ot_precond_none) THEN 00401 IF(p_env%new_preconditioner) THEN 00402 p_env%os_valid = .FALSE. 00403 DO ispin = 1,nspins 00404 CALL make_preconditioner(p_env%preconditioner(ispin),& 00405 linres_control%preconditioner_type,ot_precond_solver_default,& 00406 matrix_ks(ispin)%matrix,matrix_s(1)%matrix,matrix_t(1)%matrix,& 00407 mos(ispin)%mo_set,linres_control%energy_gap,error=error) 00408 ENDDO 00409 p_env%new_preconditioner = .FALSE. 00410 ENDIF 00411 ENDIF 00412 ! 00413 ! initalization of the linear solver 00414 ! 00415 ! A * x0 00416 CALL apply_op(qs_env,p_env,psi0_order,psi1,Ap,chc,buf,error) 00417 ! 00418 ! 00419 ! r_0 = b - Ax0 00420 DO ispin = 1,nspins 00421 CALL cp_fm_to_fm(h1_psi0(ispin)%matrix,r(ispin)%matrix,error=error) 00422 CALL cp_fm_scale_and_add(-1.0_dp,r(ispin)%matrix,-1.0_dp,Ap(ispin)%matrix,error=error) 00423 ENDDO 00424 ! 00425 ! proj r 00426 CALL postortho(r,mo_coeff_array,Sc,buf,error) 00427 ! 00428 ! preconditioner 00429 linres_control%flag="" 00430 IF(linres_control%preconditioner_type.EQ.ot_precond_none) THEN 00431 ! 00432 ! z_0 = r_0 00433 DO ispin = 1,nspins 00434 CALL cp_fm_to_fm(r(ispin)%matrix,z(ispin)%matrix,error=error) 00435 ENDDO 00436 linres_control%flag="CG" 00437 ELSE 00438 ! 00439 ! z_0 = M * r_0 00440 DO ispin = 1,nspins 00441 CALL apply_preconditioner(p_env%preconditioner(ispin),r(ispin)%matrix,& 00442 & z(ispin)%matrix,error) 00443 ENDDO 00444 linres_control%flag="PCG" 00445 ENDIF 00446 ! 00447 norm_res = 0.0_dp 00448 DO ispin = 1,nspins 00449 ! 00450 ! p_0 = z_0 00451 CALL cp_fm_to_fm(z(ispin)%matrix,p(ispin)%matrix,error=error) 00452 ! 00453 ! trace(r_0 * z_0) 00454 CALL cp_fm_trace(r(ispin)%matrix,z(ispin)%matrix,tr_rz0(ispin),error) 00455 IF(tr_rz0(ispin).LT.0.0_dp) CALL stop_program(routineN,moduleN,__LINE__,& 00456 "tr(r_j*z_j) < 0") 00457 norm_res = MAX(norm_res,ABS(tr_rz0(ispin))/SQRT(REAL(nao*maxnmo_o,dp))) 00458 ENDDO 00459 ! 00460 ! 00461 alpha(:) = 0.0_dp 00462 restart = .FALSE. 00463 should_stop = .FALSE. 00464 iteration: DO iter = 1,linres_control%max_iter 00465 ! 00466 ! check convergence 00467 linres_control%converged = .FALSE. 00468 IF(norm_res.LT.linres_control%eps) THEN 00469 linres_control%converged = .TRUE. 00470 ENDIF 00471 ! 00472 t2 = m_walltime() 00473 IF(iter.EQ.1.OR.MOD(iter,1).EQ.0.OR.linres_control%converged.OR.restart.OR.should_stop) THEN 00474 IF(output_unit>0) THEN 00475 WRITE(output_unit,"(T5,I5,T18,A3,T28,L1,T38,1E8.2,T48,F16.10,T68,F8.2)")& 00476 iter,linres_control%flag,restart,MAXVAL(alpha),norm_res,t2-t1 00477 CALL m_flush(output_unit) 00478 ENDIF 00479 ENDIF 00480 ! 00481 IF(linres_control%converged) THEN 00482 IF(output_unit>0) THEN 00483 WRITE(output_unit,"(/,T2,A,I4,A,/)") "The linear solver converged in ",iter," iterations." 00484 CALL m_flush(output_unit) 00485 ENDIF 00486 EXIT iteration 00487 ELSE IF(should_stop) THEN 00488 IF(output_unit>0) THEN 00489 WRITE(output_unit,"(/,T2,A,I4,A,/)") "The linear solver did NOT converge! External stop" 00490 CALL m_flush(output_unit) 00491 END IF 00492 EXIT iteration 00493 ENDIF 00494 ! 00495 ! Max number of iteration reached 00496 IF(iter == linres_control%max_iter) THEN 00497 IF(output_unit>0) THEN 00498 WRITE (output_unit,"(/,T2,A/)")& 00499 "The linear solver didnt converge! Maximum number of iterations reached." 00500 CALL m_flush(output_unit) 00501 ENDIF 00502 linres_control%converged = .FALSE. 00503 ENDIF 00504 ! t1 = m_walltime() 00505 ! 00506 ! 00507 ! Apply the operators that do not depend on the perturbation 00508 CALL apply_op(qs_env,p_env,psi0_order,p,Ap,chc,buf,error) 00509 ! 00510 ! 00511 ! ! proj Ap onto the virtual subspace 00512 CALL postortho(Ap,mo_coeff_array,Sc,buf,error) 00513 ! 00514 ! 00515 DO ispin = 1,nspins 00516 ! 00517 ! tr(Ap_j*p_j) 00518 CALL cp_fm_trace(Ap(ispin)%matrix,p(ispin)%matrix,tr_pAp(ispin),error) 00519 IF (tr_pAp(ispin).LT.0.0_dp) THEN 00520 CALL stop_program(routineN,moduleN,__LINE__,& 00521 "tr(Ap_j*p_j) < 0") 00522 END IF 00523 ! 00524 ! alpha = tr(r_j*z_j) / tr(Ap_j*p_j) 00525 IF(tr_pAp(ispin).LT.1.0e-10_dp) THEN 00526 alpha(ispin) = 1.0_dp 00527 ELSE 00528 alpha(ispin) = tr_rz0(ispin) / tr_pAp(ispin) 00529 ENDIF 00530 !dbg 00531 ! write(*,*) 'alpha', tr_rz0(ispin), tr_pAp(ispin), alpha(ispin) 00532 !dbg 00533 ! 00534 ! x_j+1 = x_j + alpha * p_j 00535 CALL cp_fm_scale_and_add(1.0_dp,psi1(ispin)%matrix,alpha(ispin),p(ispin)%matrix,error=error) 00536 ENDDO 00537 ! 00538 ! need to recompute the residue 00539 restart = .FALSE. 00540 IF(MOD(iter,linres_control%restart_every).EQ.0) THEN 00541 ! 00542 ! 00543 ! r_j+1 = b - A * x_j+1 00544 CALL apply_op(qs_env,p_env,psi0_order,psi1,Ap,chc,buf,error) 00545 ! 00546 ! 00547 DO ispin = 1,nspins 00548 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 00549 CALL cp_fm_to_fm(h1_psi0(ispin)%matrix,r(ispin)%matrix,error=error) 00550 CALL cp_fm_scale_and_add(-1.0_dp,r(ispin)%matrix,-1.0_dp,Ap(ispin)%matrix,error=error) 00551 ENDDO 00552 CALL postortho(r,mo_coeff_array,Sc,buf,error) 00553 ! 00554 restart = .TRUE. 00555 ELSE 00556 ! ! proj Ap onto the virtual subspace 00557 CALL postortho(Ap,mo_coeff_array,Sc,buf,error) 00558 ! 00559 ! r_j+1 = r_j - alpha * Ap_j 00560 DO ispin = 1,nspins 00561 CALL cp_fm_scale_and_add(1.0_dp,r(ispin)%matrix,-alpha(ispin),Ap(ispin)%matrix,error=error) 00562 ENDDO 00563 restart = .FALSE. 00564 ENDIF 00565 ! 00566 ! preconditioner 00567 linres_control%flag="" 00568 IF(linres_control%preconditioner_type.EQ.ot_precond_none) THEN 00569 ! 00570 ! z_j+1 = r_j+1 00571 DO ispin = 1,nspins 00572 CALL cp_fm_to_fm(r(ispin)%matrix,z(ispin)%matrix,error=error) 00573 ENDDO 00574 linres_control%flag="CG" 00575 ELSE 00576 ! 00577 ! z_j+1 = M * r_j+1 00578 DO ispin = 1,nspins 00579 CALL apply_preconditioner(p_env%preconditioner(ispin),r(ispin)%matrix,& 00580 & z(ispin)%matrix,error) 00581 ENDDO 00582 linres_control%flag="PCG" 00583 ENDIF 00584 ! 00585 norm_res = 0.0_dp 00586 DO ispin = 1,nspins 00587 ! 00588 ! tr(r_j+1*z_j+1) 00589 CALL cp_fm_trace(r(ispin)%matrix,z(ispin)%matrix,tr_rz1(ispin),error) 00590 IF(tr_rz1(ispin).LT.0.0_dp) CALL stop_program(routineN,moduleN,__LINE__,& 00591 "tr(r_j+1*z_j+1) < 0") 00592 norm_res = MAX(norm_res,tr_rz1(ispin)/SQRT(REAL(nao*maxnmo_o,dp))) 00593 ! 00594 ! beta = tr(r_j+1*z_j+1) / tr(r_j*z_j) 00595 IF(tr_rz0(ispin).LT.1.0e-10_dp) THEN 00596 beta(ispin) = 0.0_dp 00597 ELSE 00598 beta(ispin) = tr_rz1(ispin) / tr_rz0(ispin) 00599 ENDIF 00600 !dbg 00601 ! write(*,*) 'beta', tr_rz1(ispin), tr_rz0(ispin), beta(ispin) 00602 !dbg 00603 ! 00604 ! p_j+1 = z_j+1 + beta * p_j 00605 CALL cp_fm_scale_and_add(beta(ispin),p(ispin)%matrix,1.0_dp,z(ispin)%matrix,error=error) 00606 tr_rz0(ispin) = tr_rz1(ispin) 00607 ENDDO 00608 00609 ! ** can we exit the SCF loop ? 00610 CALL external_control(should_stop,"LINRES",target_time=qs_env%target_time, & 00611 start_time=qs_env%start_time,error=error) 00612 00613 00614 ENDDO iteration 00615 ! 00616 ! proj psi1 00617 CALL preortho(psi1,mo_coeff_array,Sc,buf,error) 00618 ! 00619 ! clean up 00620 DO ispin = 1,nspins 00621 CALL cp_fm_release(r(ispin)%matrix,error=error) 00622 CALL cp_fm_release(p(ispin)%matrix,error=error) 00623 CALL cp_fm_release(z(ispin)%matrix,error=error) 00624 CALL cp_fm_release(Ap(ispin)%matrix,error=error) 00625 ! 00626 CALL cp_fm_release(Sc(ispin)%matrix,error=error) 00627 CALL cp_fm_release(chc(ispin)%matrix,error=error) 00628 ENDDO 00629 CALL cp_fm_release(buf,error=error) 00630 DEALLOCATE(alpha,beta,tr_pAp,tr_rz0,tr_rz1,r,p,z,Ap,Sc,chc,mo_coeff_array,STAT=istat) 00631 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00632 ! 00633 CALL cp_print_key_finished_output(output_unit,logger,lr_section,"PRINT%PROGRAM_RUN_INFO",error=error) 00634 ! 00635 CALL timestop(handle) 00636 ! 00637 END SUBROUTINE linres_solver 00638 ! 00639 ! 00640 SUBROUTINE apply_op(qs_env,p_env,c0,v,Av,chc,buf,error) 00641 ! 00642 TYPE(qs_environment_type), POINTER :: qs_env 00643 TYPE(qs_p_env_type), POINTER :: p_env 00644 TYPE(cp_fm_p_type), DIMENSION(:), 00645 POINTER :: c0, v, Av, chc 00646 TYPE(cp_fm_type), POINTER :: buf 00647 TYPE(cp_error_type), INTENT(INOUT) :: error 00648 00649 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op', 00650 routineP = moduleN//':'//routineN 00651 00652 INTEGER :: handle, ispin, nspins 00653 LOGICAL :: failure 00654 REAL(dp) :: chksum 00655 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00656 POINTER :: matrix_ks, matrix_s 00657 TYPE(dft_control_type), POINTER :: dft_control 00658 TYPE(linres_control_type), POINTER :: linres_control 00659 TYPE(qs_rho_type), POINTER :: rho 00660 00661 ! 00662 00663 failure = .FALSE. 00664 NULLIFY(dft_control,matrix_ks,matrix_s,linres_control) 00665 ! 00666 ! 00667 CALL timeset(routineN,handle) 00668 ! 00669 CPPrecondition(ASSOCIATED(v),cp_failure_level,routineP,error,failure) 00670 CPPrecondition(ASSOCIATED(Av),cp_failure_level,routineP,error,failure) 00671 CPPrecondition(ASSOCIATED(chc),cp_failure_level,routineP,error,failure) 00672 CPPrecondition(ASSOCIATED(buf),cp_failure_level,routineP,error,failure) 00673 ! 00674 CALL get_qs_env(qs_env=qs_env,& 00675 matrix_ks=matrix_ks,& 00676 matrix_s=matrix_s,& 00677 dft_control=dft_control,& 00678 linres_control=linres_control,& 00679 error=error) 00680 00681 ! 00682 nspins = dft_control%nspins 00683 ! 00684 ! 00685 ! apply the uncoupled operator 00686 DO ispin = 1,nspins 00687 CALL apply_op_1(v(ispin)%matrix,Av(ispin)%matrix,matrix_ks(ispin)%matrix,& 00688 matrix_s(1)%matrix,chc(ispin)%matrix,buf,error) 00689 ENDDO 00690 00691 00692 IF(linres_control%do_kernel) THEN!.AND.chksum.GT.1.0E-10_dp) THEN 00693 ! 00694 ! build DM 00695 CALL build_dm_response(c0,v,p_env%p1,error) 00696 00697 chksum = 0.0_dp 00698 DO ispin = 1,nspins 00699 chksum = chksum + cp_dbcsr_checksum(p_env%p1(ispin)%matrix,error=error) 00700 ENDDO 00701 00702 ! 00703 ! skip the kernel if the DM is very small 00704 IF(chksum.GT.1.0E-14_dp) THEN 00705 00706 CALL p_env_check_i_alloc(p_env,qs_env,error) 00707 00708 DO ispin = 1,nspins 00709 CALL cp_dbcsr_copy(p_env%rho1%rho_ao(ispin)%matrix,p_env%p1(ispin)%matrix,error=error) 00710 ENDDO 00711 00712 CALL qs_rho_update_rho(rho_struct=p_env%rho1,local_rho_set=p_env%local_rho_set,& 00713 qs_env=qs_env,error=error) 00714 00715 !if(first_time) then 00716 CALL get_qs_env(qs_env,rho=rho,error=error) ! that could be called before 00717 CALL qs_rho_update_rho(rho,qs_env=qs_env,error=error) ! that could be called before 00718 ! first_time = .false. 00719 !endif 00720 00721 00722 00723 CALL apply_op_2(qs_env,p_env,c0,v,Av,chc,buf,error) 00724 00725 !CALL kpp1_calc_k_p_p1(p_env%kpp1_env, p_env, qs_env, p_env%kpp1, qs_env%rho, p_env%rho1, p_env%rho1_xc, error) 00726 !DO ispin = 1,nspins 00727 ! CALL cp_fm_get_info(c0(ispin)%matrix,ncol_global=ncol,error=error) 00728 ! CALL cp_sm_fm_multiply(sparse_matrix=p_env%kpp1(ispin)%matrix,& 00729 ! v_in=c0(ispin)%matrix,& 00730 ! v_out=Av(ispin)%matrix,& 00731 ! ncol=ncol,alpha=1.0_dp,beta=1.0_dp,& 00732 ! error=error) 00733 !ENDDO 00734 00735 ENDIF 00736 00737 ENDIF 00738 ! 00739 CALL timestop(handle) 00740 ! 00741 END SUBROUTINE apply_op 00742 ! 00743 ! 00744 SUBROUTINE build_dm_response(c0,c1,dm,error) 00745 ! 00746 TYPE(cp_fm_p_type), DIMENSION(:), 00747 POINTER :: c0, c1 00748 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00749 POINTER :: dm 00750 TYPE(cp_error_type), INTENT(INOUT) :: error 00751 00752 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dm_response', 00753 routineP = moduleN//':'//routineN 00754 00755 INTEGER :: ispin, ncol, nspins 00756 LOGICAL :: failure 00757 00758 ! 00759 ! 00760 00761 failure = .FALSE. 00762 00763 CPPrecondition(ASSOCIATED(c0),cp_failure_level,routineP,error,failure) 00764 CPPrecondition(ASSOCIATED(c1),cp_failure_level,routineP,error,failure) 00765 CPPrecondition(ASSOCIATED(dm),cp_failure_level,routineP,error,failure) 00766 00767 nspins = SIZE(dm,1) 00768 00769 DO ispin = 1,nspins 00770 CALL cp_dbcsr_set(dm(ispin)%matrix,0.0_dp,error=error) 00771 CALL cp_fm_get_info(c0(ispin)%matrix,ncol_global=ncol,error=error) 00772 CALL cp_dbcsr_plus_fm_fm_t(dm(ispin)%matrix,& 00773 matrix_v=c0(ispin)%matrix,& 00774 matrix_g=c1(ispin)%matrix,& 00775 ncol=ncol,alpha=1.0_dp,error=error) 00776 CALL cp_dbcsr_plus_fm_fm_t(dm(ispin)%matrix,& 00777 matrix_v=c1(ispin)%matrix,& 00778 matrix_g=c0(ispin)%matrix,& 00779 ncol=ncol,alpha=1.0_dp,error=error) 00780 ENDDO 00781 00782 END SUBROUTINE build_dm_response 00783 ! 00784 ! 00785 SUBROUTINE apply_op_1(v,Av,matrix_ks,matrix_s,chc,buf,error) 00786 ! 00787 TYPE(cp_fm_type), POINTER :: v, Av 00788 TYPE(cp_dbcsr_type), POINTER :: matrix_ks, matrix_s 00789 TYPE(cp_fm_type), POINTER :: chc, buf 00790 TYPE(cp_error_type), INTENT(INOUT) :: error 00791 00792 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op_1', 00793 routineP = moduleN//':'//routineN 00794 00795 INTEGER :: handle, ncol, nrow 00796 LOGICAL :: failure 00797 00798 ! 00799 00800 failure = .FALSE. 00801 ! 00802 CALL timeset(routineN,handle) 00803 ! 00804 CPPrecondition(ASSOCIATED(v),cp_failure_level,routineP,error,failure) 00805 CPPrecondition(ASSOCIATED(Av),cp_failure_level,routineP,error,failure) 00806 CPPrecondition(ASSOCIATED(matrix_ks),cp_failure_level,routineP,error,failure) 00807 CPPrecondition(ASSOCIATED(matrix_s),cp_failure_level,routineP,error,failure) 00808 CPPrecondition(ASSOCIATED(chc),cp_failure_level,routineP,error,failure) 00809 CPPrecondition(ASSOCIATED(buf),cp_failure_level,routineP,error,failure) 00810 ! 00811 CALL cp_fm_get_info(v,ncol_global=ncol,nrow_global=nrow,error=error) 00812 ! H * v 00813 CALL cp_dbcsr_sm_fm_multiply(matrix_ks,v,Av,ncol,error=error) 00814 ! v * e (chc already multiplied by -1) 00815 CALL cp_fm_gemm('N','N',nrow,ncol,ncol,1.0_dp,v,chc,0.0_dp,buf,error) 00816 ! S * ve 00817 CALL cp_dbcsr_sm_fm_multiply(matrix_s,buf,Av,ncol,alpha=1.0_dp,beta=1.0_dp,error=error) 00818 !Results is H*C1 - S*<iHj>*C1 00819 ! 00820 CALL timestop(handle) 00821 ! 00822 END SUBROUTINE apply_op_1 00823 00824 00825 SUBROUTINE apply_op_2(qs_env,p_env,c0,v,Av,chc,buf,error) 00826 ! 00827 TYPE(qs_environment_type), POINTER :: qs_env 00828 TYPE(qs_p_env_type), POINTER :: p_env 00829 TYPE(cp_fm_p_type), DIMENSION(:), 00830 POINTER :: c0, v, Av, chc 00831 TYPE(cp_fm_type), POINTER :: buf 00832 TYPE(cp_error_type), INTENT(INOUT) :: error 00833 00834 CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2', 00835 routineP = moduleN//':'//routineN 00836 00837 INTEGER :: handle, ispin, ncol, nspins, 00838 stat 00839 INTEGER, DIMENSION(2, 3) :: bo 00840 LOGICAL :: failure, gapw, gapw_xc, 00841 lr_triplet, lsd 00842 REAL(KIND=dp) :: energy_hartree, 00843 energy_hartree_1c, fac 00844 TYPE(cell_type), POINTER :: cell 00845 TYPE(cp_logger_type), POINTER :: logger 00846 TYPE(dft_control_type), POINTER :: dft_control 00847 TYPE(linres_control_type), POINTER :: linres_control 00848 TYPE(pw_env_type), POINTER :: pw_env 00849 TYPE(pw_p_type) :: rho1_tot_gspace, 00850 v_hartree_gspace, 00851 v_hartree_rspace 00852 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g_pw, rho1_r, rho1_r_pw, 00853 tau_pw, v_rspace_new, v_xc 00854 TYPE(pw_poisson_type), POINTER :: poisson_env 00855 TYPE(pw_pool_p_type), DIMENSION(:), 00856 POINTER :: pw_pools 00857 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00858 TYPE(qs_rho_type), POINTER :: rho, rho1, rho1_xc 00859 TYPE(section_vals_type), POINTER :: input, scf_section, 00860 xc_fun_section, xc_section 00861 TYPE(xc_rho_cflags_type) :: needs 00862 TYPE(xc_rho_set_type), POINTER :: rho1_set 00863 00864 CALL timeset(routineN,handle) 00865 failure=.FALSE. 00866 00867 NULLIFY(auxbas_pw_pool, pw_pools, pw_env, v_rspace_new, & 00868 rho1_r, rho1_g_pw, tau_pw, cell, v_xc, rho1_set,& 00869 poisson_env, input, scf_section,rho,dft_control,logger) 00870 logger => cp_error_get_logger(error) 00871 00872 energy_hartree=0.0_dp 00873 energy_hartree_1c=0.0_dp 00874 00875 CPPrecondition(ASSOCIATED(c0),cp_failure_level,routineP,error,failure) 00876 CPPrecondition(ASSOCIATED(v),cp_failure_level,routineP,error,failure) 00877 CPPrecondition(ASSOCIATED(Av),cp_failure_level,routineP,error,failure) 00878 CPPrecondition(ASSOCIATED(chc),cp_failure_level,routineP,error,failure) 00879 CPPrecondition(ASSOCIATED(buf),cp_failure_level,routineP,error,failure) 00880 00881 00882 CPPrecondition(ASSOCIATED(p_env%kpp1_env),cp_failure_level,routineP,error,failure) 00883 CPPrecondition(ASSOCIATED(p_env%kpp1),cp_failure_level,routineP,error,failure) 00884 rho1 => p_env%rho1 00885 rho1_xc => p_env%rho1_xc 00886 CPPrecondition(ASSOCIATED(rho1),cp_failure_level,routineP,error,failure) 00887 00888 IF (.NOT.failure) THEN 00889 CPPrecondition(p_env%kpp1_env%ref_count>0,cp_failure_level,routineP,error,failure) 00890 00891 CALL get_qs_env(qs_env=qs_env,& 00892 pw_env=pw_env,& 00893 cell=cell,& 00894 input=input,& 00895 rho=rho,& 00896 linres_control=linres_control,& 00897 dft_control=dft_control,& 00898 error=error) 00899 00900 lr_triplet = linres_control%lr_triplet 00901 CALL kpp1_check_i_alloc(p_env%kpp1_env,qs_env,lr_triplet,error=error) 00902 !gapw=(section_get_ival(input,"DFT%QS%METHOD",error=error)==do_method_gapw) 00903 !gapw_xc=(section_get_ival(input,"DFT%QS%METHOD",error=error)==do_method_gapw_xc) 00904 gapw = dft_control%qs_control%gapw 00905 gapw_xc = dft_control%qs_control%gapw_xc 00906 IF(gapw_xc) THEN 00907 CPPrecondition(ASSOCIATED(rho1_xc),cp_failure_level,routineP,error,failure) 00908 END IF 00909 00910 nspins = SIZE(p_env%kpp1) 00911 lsd = (nspins==2) 00912 00913 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 00914 scf_section => section_vals_get_subs_vals(input,"DFT%SCF",error=error) 00915 END IF 00916 00917 IF (.NOT.failure) THEN 00918 p_env%kpp1_env%iter=p_env%kpp1_env%iter+1 00919 END IF 00920 00921 ! gets the tmp grids 00922 IF (.NOT. failure) THEN 00923 CPPrecondition(ASSOCIATED(pw_env),cp_failure_level,routineP,error,failure) 00924 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 00925 pw_pools=pw_pools, poisson_env=poisson_env,error=error) 00926 ALLOCATE(v_rspace_new(nspins), stat=stat) 00927 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00928 END IF 00929 IF (.NOT.failure) THEN 00930 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_gspace%pw,& 00931 use_data = COMPLEXDATA1D,& 00932 in_space = RECIPROCALSPACE, error=error) 00933 CALL pw_pool_create_pw(auxbas_pw_pool,v_hartree_rspace%pw,& 00934 use_data = REALDATA3D,& 00935 in_space = REALSPACE, error=error) 00936 END IF 00937 00938 IF (gapw .OR. gapw_xc) & 00939 CALL prepare_gapw_den(qs_env,p_env%local_rho_set, do_rho0=(.NOT.gapw_xc), error=error) 00940 00941 ! *** calculate the hartree potential on the total density *** 00942 IF (.NOT. failure) THEN 00943 00944 CALL pw_pool_create_pw(auxbas_pw_pool, rho1_tot_gspace%pw,& 00945 use_data = COMPLEXDATA1D,& 00946 in_space = RECIPROCALSPACE, error=error) 00947 00948 CALL pw_copy(rho1%rho_g(1)%pw,rho1_tot_gspace%pw, error=error) 00949 DO ispin=2,nspins 00950 CALL pw_axpy(rho1%rho_g(ispin)%pw, rho1_tot_gspace%pw, error=error) 00951 END DO 00952 IF (gapw) & 00953 CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs%pw, rho1_tot_gspace%pw,& 00954 error=error) 00955 00956 !IF (cp_print_key_should_output(logger%iter_info,scf_section,"PRINT%TOTAL_DENSITIES",& 00957 ! error=error)/=0) THEN 00958 ! output_unit = cp_print_key_unit_nr(logger,scf_section,"PRINT%TOTAL_DENSITIES",& 00959 ! extension=".scfLog",error=error) 00960 ! CALL print_densities(kpp1_env, rho1, rho1_tot_gspace, output_unit, error=error) 00961 ! CALL cp_print_key_finished_output(output_unit,logger,scf_section,& 00962 ! "PRINT%TOTAL_DENSITIES", error=error) 00963 !END IF 00964 00965 IF (.NOT.(nspins==1 .AND. lr_triplet )) THEN 00966 CALL pw_poisson_solve(poisson_env,rho1_tot_gspace%pw, & 00967 energy_hartree, & 00968 v_hartree_gspace%pw,error=error) 00969 CALL pw_transfer(v_hartree_gspace%pw,v_hartree_rspace%pw, error=error) 00970 ENDIF 00971 00972 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho1_tot_gspace%pw,& 00973 error=error) 00974 00975 ! *** calculate the xc potential *** 00976 IF(gapw_xc) THEN 00977 CALL qs_rho_get(rho1_xc, rho_r=rho1_r, error=error) 00978 ELSE 00979 CALL qs_rho_get(rho1, rho_r=rho1_r, error=error) 00980 END IF 00981 00982 IF (nspins == 1 .AND. (lr_triplet)) THEN 00983 00984 lsd = .TRUE. 00985 ALLOCATE(rho1_r_pw(2)) 00986 DO ispin=1, 2 00987 NULLIFY(rho1_r_pw(ispin)%pw) 00988 CALL pw_create(rho1_r_pw(ispin)%pw, rho1_r(1)%pw%pw_grid, & 00989 rho1_r(1)%pw%in_use, rho1_r(1)%pw%in_space,error=error) 00990 CALL pw_transfer(rho1_r(1)%pw, rho1_r_pw(ispin)%pw, error=error) 00991 END DO 00992 00993 ELSE 00994 00995 ALLOCATE(rho1_r_pw(nspins)) 00996 DO ispin=1, nspins 00997 rho1_r_pw(ispin)%pw => rho1_r(ispin)%pw 00998 CALL pw_retain(rho1_r_pw(ispin)%pw,error=error) 00999 END DO 01000 01001 END IF 01002 01003 NULLIFY(tau_pw) 01004 01005 !------! 01006 ! rho1 ! 01007 !------! 01008 bo = rho1_r(1)%pw%pw_grid%bounds_local 01009 ! create the place where to store the argument for the functionals 01010 CALL xc_rho_set_create(rho1_set, bo, & 01011 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF",error=error), & 01012 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF",error=error), & 01013 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF",error=error), & 01014 error=error) 01015 01016 xc_fun_section => section_vals_get_subs_vals(xc_section,"XC_FUNCTIONAL",& 01017 error=error) 01018 needs=xc_functionals_get_needs(xc_fun_section,lsd,.TRUE.,error) 01019 01020 ! calculate the arguments needed by the functionals 01021 CALL xc_rho_set_update(rho1_set, rho1_r_pw, rho1_g_pw, tau_pw, needs,& 01022 section_get_ival(xc_section,"XC_GRID%XC_DERIV",error=error),& 01023 section_get_ival(xc_section,"XC_GRID%XC_SMOOTH_RHO",error=error),& 01024 cell, auxbas_pw_pool, error) 01025 01026 ALLOCATE(v_xc(nspins),stat=stat) 01027 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01028 DO ispin=1, nspins 01029 NULLIFY(v_xc(ispin)%pw) 01030 CALL pw_pool_create_pw(auxbas_pw_pool,v_xc(ispin)%pw,& 01031 use_data = REALDATA3D,& 01032 in_space = REALSPACE, error=error) 01033 CALL pw_zero(v_xc(ispin)%pw, error=error) 01034 END DO 01035 01036 fac=0._dp 01037 IF (nspins==1) THEN 01038 IF (lr_triplet) fac=-1.0_dp 01039 END IF 01040 01041 CALL xc_calc_2nd_deriv(v_xc, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, & 01042 rho1_set, auxbas_pw_pool,xc_section=xc_section,& 01043 tddfpt_fac=fac, error=error) 01044 01045 DO ispin=1,nspins 01046 v_rspace_new(ispin)%pw => v_xc(ispin)%pw 01047 END DO 01048 DEALLOCATE(v_xc,stat=stat) 01049 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01050 01051 IF(gapw) CALL calculate_xc_2nd_deriv_atom(p_env,qs_env,xc_section,do_tddft=.FALSE.,& 01052 do_triplet=lr_triplet,error=error) 01053 01054 CALL xc_rho_set_release(rho1_set,error=error) 01055 01056 DO ispin=1,SIZE(rho1_r_pw) 01057 CALL pw_release(rho1_r_pw(ispin)%pw,error=error) 01058 END DO 01059 DEALLOCATE(rho1_r_pw, stat=stat) 01060 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01061 01062 !-------------------------------! 01063 ! Add both hartree and xc terms ! 01064 !-------------------------------! 01065 DO ispin=1,nspins 01066 01067 IF(gapw_xc) THEN 01068 ! XC and Hartree are integrated separatedly 01069 ! XC uses the sofft basis set only 01070 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d * & 01071 v_rspace_new(ispin)%pw%pw_grid%dvol 01072 01073 IF ( nspins==1) THEN 01074 01075 IF (.NOT.(lr_triplet)) THEN 01076 01077 v_rspace_new(1)%pw%cr3d = 2.0_dp * v_rspace_new(1)%pw%cr3d 01078 01079 END IF 01080 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 01081 CALL cp_dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix,0.0_dp,error=error) 01082 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 01083 p=rho%rho_ao(ispin),& 01084 h=p_env%kpp1_env%v_ao(ispin),& 01085 qs_env=qs_env,& 01086 calculate_forces=.FALSE.,gapw=gapw_xc,error=error) 01087 01088 ! add hartree only for SINGLETS 01089 IF (.NOT.lr_triplet) THEN 01090 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 01091 v_hartree_rspace%pw%pw_grid%dvol 01092 v_rspace_new(1)%pw%cr3d = 2.0_dp * v_hartree_rspace%pw%cr3d 01093 01094 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 01095 p=rho%rho_ao(ispin),& 01096 h=p_env%kpp1_env%v_ao(ispin),& 01097 qs_env=qs_env,& 01098 calculate_forces=.FALSE.,gapw=gapw,error=error) 01099 END IF 01100 ELSE 01101 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 01102 CALL cp_dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix,0.0_dp,error=error) 01103 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 01104 p=rho%rho_ao(ispin),& 01105 h=p_env%kpp1_env%v_ao(ispin),& 01106 qs_env=qs_env,& 01107 calculate_forces=.FALSE.,gapw=gapw_xc,error=error) 01108 01109 IF (ispin == 1) THEN 01110 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 01111 v_hartree_rspace%pw%pw_grid%dvol 01112 END IF 01113 v_rspace_new(ispin)%pw%cr3d = v_hartree_rspace%pw%cr3d 01114 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 01115 p=rho%rho_ao(ispin),& 01116 h=p_env%kpp1_env%v_ao(ispin),& 01117 qs_env=qs_env,& 01118 calculate_forces=.FALSE.,gapw=gapw,error=error) 01119 END IF 01120 01121 ELSE 01122 01123 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d * & 01124 v_rspace_new(ispin)%pw%pw_grid%dvol 01125 01126 IF ( nspins==1) THEN 01127 01128 IF(.NOT.(lr_triplet)) THEN 01129 v_rspace_new(1)%pw%cr3d = 2.0_dp * v_rspace_new(1)%pw%cr3d 01130 01131 ENDIF 01132 01133 ! add hartree only for SINGLETS 01134 !IF (res_etype == tddfpt_singlet) THEN 01135 IF (.NOT.lr_triplet) THEN 01136 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 01137 v_hartree_rspace%pw%pw_grid%dvol 01138 v_rspace_new(1)%pw%cr3d = v_rspace_new(1)%pw%cr3d + & 01139 2.0_dp * v_hartree_rspace%pw%cr3d 01140 END IF 01141 ELSE 01142 IF (ispin == 1) THEN 01143 v_hartree_rspace%pw%cr3d = v_hartree_rspace%pw%cr3d * & 01144 v_hartree_rspace%pw%pw_grid%dvol 01145 END IF 01146 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d + & 01147 v_hartree_rspace%pw%cr3d 01148 END IF 01149 01150 ! remove kpp1_env%v_ao and work directly on k_p_p1 ? 01151 CALL cp_dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix,0.0_dp,error=error) 01152 01153 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin),& 01154 p=rho%rho_ao(ispin),& 01155 h=p_env%kpp1_env%v_ao(ispin),& 01156 qs_env=qs_env,& 01157 calculate_forces=.FALSE.,gapw=gapw,error=error) 01158 01159 END IF 01160 01161 CALL cp_dbcsr_copy(p_env%kpp1(ispin)%matrix,p_env%kpp1_env%v_ao(ispin)%matrix,& 01162 error=error) 01163 END DO 01164 01165 IF (gapw) THEN 01166 01167 IF (.NOT. ( (nspins == 1 .AND. lr_triplet))) THEN 01168 CALL Vh_1c_gg_integrals(qs_env,energy_hartree_1c, tddft=.TRUE., do_triplet=lr_triplet, & 01169 p_env=p_env,error=error) 01170 01171 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, & 01172 .FALSE., tddft=.TRUE., do_triplet=lr_triplet, & 01173 p_env=p_env, error=error) 01174 END IF 01175 01176 ! *** Add single atom contributions to the KS matrix *** 01177 01178 CALL update_ks_atom(qs_env,p_env%kpp1,rho%rho_ao,.FALSE.,.TRUE.,p_env,error) 01179 01180 END IF 01181 01182 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_gspace%pw,& 01183 error=error) 01184 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_hartree_rspace%pw,& 01185 error=error) 01186 DO ispin=1,nspins 01187 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_rspace_new(ispin)%pw,& 01188 error=error) 01189 END DO 01190 DEALLOCATE(v_rspace_new, stat=stat) 01191 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01192 01193 END IF 01194 01195 DO ispin = 1,nspins 01196 CALL cp_fm_get_info(c0(ispin)%matrix,ncol_global=ncol,error=error) 01197 CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1(ispin)%matrix,& 01198 c0(ispin)%matrix,& 01199 Av(ispin)%matrix,& 01200 ncol=ncol,alpha=1.0_dp,beta=1.0_dp,& 01201 error=error) 01202 ENDDO 01203 ! 01204 CALL timestop(handle) 01205 ! 01206 END SUBROUTINE apply_op_2 01207 ! 01208 ! 01209 SUBROUTINE p_env_check_i_alloc(p_env, qs_env, error) 01210 TYPE(qs_p_env_type), POINTER :: p_env 01211 TYPE(qs_environment_type), POINTER :: qs_env 01212 TYPE(cp_error_type), INTENT(inout) :: error 01213 01214 CHARACTER(len=*), PARAMETER :: routineN = 'p_env_check_i_alloc', 01215 routineP = moduleN//':'//routineN 01216 01217 CHARACTER(len=25) :: name 01218 INTEGER :: handle, ispin, nspins 01219 LOGICAL :: failure, gapw_xc 01220 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01221 POINTER :: matrix_s 01222 TYPE(dft_control_type), POINTER :: dft_control 01223 01224 CALL timeset(routineN,handle) 01225 01226 failure=.FALSE. 01227 NULLIFY(dft_control,matrix_s) 01228 01229 CPPrecondition(ASSOCIATED(p_env),cp_failure_level,routineP,error,failure) 01230 CPPrecondition(p_env%ref_count>0,cp_failure_level,routineP,error,failure) 01231 IF (.NOT. failure) THEN 01232 CALL get_qs_env(qs_env, dft_control=dft_control,error=error) 01233 gapw_xc = dft_control%qs_control%gapw_xc 01234 IF (.NOT.ASSOCIATED(p_env%kpp1)) THEN 01235 CALL get_qs_env(qs_env, matrix_s=matrix_s, error=error) 01236 nspins=dft_control%nspins 01237 01238 CALL cp_dbcsr_allocate_matrix_set(p_env%kpp1,nspins,error=error) 01239 name="p_env"//cp_to_string(p_env%id_nr)//"%kpp1-" 01240 !CALL compress(name,full=.TRUE.) 01241 DO ispin=1,nspins 01242 ALLOCATE(p_env%kpp1(ispin)%matrix) 01243 CALL cp_dbcsr_init(p_env%kpp1(ispin)%matrix,error=error) 01244 CALL cp_dbcsr_copy(p_env%kpp1(ispin)%matrix,matrix_s(1)%matrix,& 01245 name=TRIM(name)//ADJUSTL(cp_to_string(ispin)),error=error) 01246 CALL cp_dbcsr_set(p_env%kpp1(ispin)%matrix,0.0_dp,error=error) 01247 END DO 01248 01249 CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env, error=error) 01250 IF(gapw_xc) THEN 01251 CALL qs_rho_rebuild(p_env%rho1_xc,qs_env=qs_env,gapw_xc=gapw_xc,error=error) 01252 p_env%rho1_xc%rho_ao => p_env%rho1%rho_ao 01253 END IF 01254 END IF 01255 01256 IF (.NOT.ASSOCIATED(p_env%rho1)) THEN 01257 CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env, error=error) 01258 IF(gapw_xc) THEN 01259 CALL qs_rho_rebuild(p_env%rho1_xc,qs_env=qs_env,gapw_xc=gapw_xc,error=error) 01260 p_env%rho1_xc%rho_ao => p_env%rho1%rho_ao 01261 END IF 01262 END IF 01263 END IF 01264 CALL timestop(handle) 01265 END SUBROUTINE p_env_check_i_alloc 01266 ! 01267 ! 01268 SUBROUTINE kpp1_check_i_alloc(kpp1_env, qs_env, lr_triplet,error) 01269 01270 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env 01271 TYPE(qs_environment_type), POINTER :: qs_env 01272 LOGICAL, INTENT(IN) :: lr_triplet 01273 TYPE(cp_error_type), INTENT(inout) :: error 01274 01275 CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_check_i_alloc', 01276 routineP = moduleN//':'//routineN 01277 01278 INTEGER :: ispin, nspins, stat 01279 LOGICAL :: failure 01280 TYPE(cell_type), POINTER :: cell 01281 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01282 POINTER :: matrix_s 01283 TYPE(pw_env_type), POINTER :: pw_env 01284 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r 01285 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01286 TYPE(qs_rho_type), POINTER :: rho 01287 TYPE(section_vals_type), POINTER :: input, xc_section 01288 01289 ! ------------------------------------------------------------------ 01290 01291 failure=.FALSE. 01292 NULLIFY(pw_env,auxbas_pw_pool,matrix_s,rho,rho_r,input) 01293 01294 CPPrecondition(ASSOCIATED(kpp1_env),cp_failure_level,routineP,error,failure) 01295 CPPrecondition(kpp1_env%ref_count>0,cp_failure_level,routineP,error,failure) 01296 01297 CALL get_qs_env(qs_env, pw_env=pw_env,& 01298 matrix_s=matrix_s, cell=cell, input=input, error=error, rho=rho) 01299 01300 CALL qs_rho_get(rho, rho_r=rho_r, error=error) 01301 nspins=SIZE(rho_r) 01302 01303 CALL pw_env_get(pw_env, auxbas_pw_pool = auxbas_pw_pool, error=error) 01304 01305 IF (.NOT.ASSOCIATED(kpp1_env%v_rspace)) THEN 01306 ALLOCATE(kpp1_env%v_rspace(nspins),stat=stat) 01307 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01308 IF (.NOT.failure) THEN 01309 DO ispin=1,nspins 01310 CALL pw_pool_create_pw(auxbas_pw_pool, & 01311 kpp1_env%v_rspace(ispin)%pw,& 01312 use_data=REALDATA3D, in_space=REALSPACE,error=error) 01313 END DO 01314 END IF 01315 END IF 01316 01317 IF (.NOT.ASSOCIATED(kpp1_env%v_ao)) THEN 01318 CALL cp_dbcsr_allocate_matrix_set(kpp1_env%v_ao,nspins,error) 01319 DO ispin=1,nspins 01320 ALLOCATE(kpp1_env%v_ao(ispin)%matrix) 01321 CALL cp_dbcsr_init(kpp1_env%v_ao(ispin)%matrix,error=error) 01322 CALL cp_dbcsr_copy(kpp1_env%v_ao(ispin)%matrix,matrix_s(1)%matrix,& 01323 name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)),error=error) 01324 END DO 01325 END IF 01326 01327 IF(.NOT.ASSOCIATED(kpp1_env%deriv_set)) THEN 01328 01329 IF (nspins==1.AND.lr_triplet) THEN 01330 ALLOCATE(my_rho_r(2),stat=stat) 01331 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01332 DO ispin=1,2 01333 CALL pw_pool_create_pw(auxbas_pw_pool,my_rho_r(ispin)%pw, & 01334 use_data=rho_r(1)%pw%in_use, in_space=rho_r(1)%pw%in_space,& 01335 error=error) 01336 my_rho_r(ispin)%pw%cr3d = 0.5_dp * rho_r(1)%pw%cr3d 01337 END DO 01338 ELSE 01339 ALLOCATE(my_rho_r(SIZE(rho_r)),stat=stat) 01340 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01341 DO ispin=1,SIZE(rho_r) 01342 my_rho_r(ispin)%pw => rho_r(ispin)%pw 01343 CALL pw_retain(my_rho_r(ispin)%pw,error=error) 01344 END DO 01345 END IF 01346 01347 !ALLOCATE(my_rho_r(SIZE(rho_r)),stat=stat) 01348 !CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01349 !DO ispin=1,SIZE(rho_r) 01350 ! my_rho_r(ispin)%pw => rho_r(ispin)%pw 01351 ! CALL pw_retain(my_rho_r(ispin)%pw,error=error) 01352 !END DO 01353 01354 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 01355 01356 CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, & 01357 my_rho_r, auxbas_pw_pool, & 01358 xc_section=xc_section, cell=cell, error=error) 01359 01360 DO ispin=1,SIZE(my_rho_r) 01361 CALL pw_release(my_rho_r(ispin)%pw,error=error) 01362 ENDDO 01363 DEALLOCATE(my_rho_r,stat=stat) 01364 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01365 ENDIF 01366 01367 END SUBROUTINE kpp1_check_i_alloc 01368 ! 01369 ! 01370 SUBROUTINE preortho(v,psi0,S_psi0,buf,error) 01371 !v = (I-PS)v 01372 ! 01373 TYPE(cp_fm_p_type), DIMENSION(:), 01374 POINTER :: v, psi0, S_psi0 01375 TYPE(cp_fm_type), POINTER :: buf 01376 TYPE(cp_error_type), INTENT(INOUT) :: error 01377 01378 CHARACTER(LEN=*), PARAMETER :: routineN = 'preortho', 01379 routineP = moduleN//':'//routineN 01380 01381 INTEGER :: handle, ispin, mp, mt, mv, 01382 np, nspins, nt, nv 01383 LOGICAL :: failure 01384 01385 ! 01386 01387 failure = .FALSE. 01388 ! 01389 CALL timeset(routineN,handle) 01390 ! 01391 CPPrecondition(ASSOCIATED(v),cp_failure_level,routineP,error,failure) 01392 CPPrecondition(ASSOCIATED(S_psi0),cp_failure_level,routineP,error,failure) 01393 CPPrecondition(ASSOCIATED(psi0),cp_failure_level,routineP,error,failure) 01394 CPPrecondition(ASSOCIATED(buf),cp_failure_level,routineP,error,failure) 01395 ! 01396 nspins = SIZE(v,1) 01397 ! 01398 DO ispin = 1,nspins 01399 ! 01400 CALL cp_fm_get_info(v(ispin)%matrix,ncol_global=mv,nrow_global=nv,error=error) 01401 CALL cp_fm_get_info(psi0(ispin)%matrix,ncol_global=mp,nrow_global=np,error=error) 01402 CALL cp_fm_get_info(buf,ncol_global=mt,nrow_global=nt,error=error) 01403 CPPrecondition(nv==np,cp_failure_level,routineP,error,failure) 01404 CPPrecondition(mt>=mv,cp_failure_level,routineP,error,failure) 01405 CPPrecondition(mt>=mp,cp_failure_level,routineP,error,failure) 01406 CPPrecondition(nt==nv,cp_failure_level,routineP,error,failure) 01407 ! 01408 ! buf = v' * S_psi0 01409 CALL cp_fm_gemm('T','N',mv,mp,nv,1.0_dp,v(ispin)%matrix,S_psi0(ispin)%matrix,0.0_dp,buf,error) 01410 ! v = v - psi0 * buf' 01411 CALL cp_fm_gemm('N','T',nv,mv,mp,-1.0_dp,psi0(ispin)%matrix,buf,1.0_dp,v(ispin)%matrix,error) 01412 ! 01413 ENDDO 01414 ! 01415 CALL timestop(handle) 01416 ! 01417 END SUBROUTINE preortho 01418 ! 01419 ! 01420 SUBROUTINE postortho(v,psi0,S_psi0,buf,error) 01421 !v = (I-SP)v 01422 ! 01423 TYPE(cp_fm_p_type), DIMENSION(:), 01424 POINTER :: v, psi0, S_psi0 01425 TYPE(cp_fm_type), POINTER :: buf 01426 TYPE(cp_error_type), INTENT(INOUT) :: error 01427 01428 CHARACTER(LEN=*), PARAMETER :: routineN = 'postortho', 01429 routineP = moduleN//':'//routineN 01430 01431 INTEGER :: handle, ispin, mp, mt, mv, 01432 np, nspins, nt, nv 01433 LOGICAL :: failure 01434 01435 failure = .FALSE. 01436 ! 01437 CALL timeset(routineN,handle) 01438 ! 01439 CPPrecondition(ASSOCIATED(v),cp_failure_level,routineP,error,failure) 01440 CPPrecondition(ASSOCIATED(S_psi0),cp_failure_level,routineP,error,failure) 01441 CPPrecondition(ASSOCIATED(psi0),cp_failure_level,routineP,error,failure) 01442 CPPrecondition(ASSOCIATED(buf),cp_failure_level,routineP,error,failure) 01443 ! 01444 nspins = SIZE(v,1) 01445 ! 01446 DO ispin = 1,nspins 01447 ! 01448 CALL cp_fm_get_info(v(ispin)%matrix,ncol_global=mv,nrow_global=nv,error=error) 01449 CALL cp_fm_get_info(psi0(ispin)%matrix,ncol_global=mp,nrow_global=np,error=error) 01450 CALL cp_fm_get_info(buf,ncol_global=mt,nrow_global=nt,error=error) 01451 CPPrecondition(nv==np,cp_failure_level,routineP,error,failure) 01452 CPPrecondition(mt>=mv,cp_failure_level,routineP,error,failure) 01453 CPPrecondition(mt>=mp,cp_failure_level,routineP,error,failure) 01454 CPPrecondition(nt==nv,cp_failure_level,routineP,error,failure) 01455 ! 01456 ! buf = v' * psi0 01457 CALL cp_fm_gemm('T','N',mv,mp,nv,1.0_dp,v(ispin)%matrix,psi0(ispin)%matrix,0.0_dp,buf,error) 01458 ! v = v - S_psi0 * buf' 01459 CALL cp_fm_gemm('N','T',nv,mv,mp,-1.0_dp,S_psi0(ispin)%matrix,buf,1.0_dp,v(ispin)%matrix,error) 01460 ! 01461 ENDDO 01462 ! 01463 CALL timestop(handle) 01464 ! 01465 END SUBROUTINE postortho 01466 01467 SUBROUTINE linres_write_restart(qs_env,linres_section,vec,ivec,tag,ind,error) 01468 TYPE(qs_environment_type), POINTER :: qs_env 01469 TYPE(section_vals_type), POINTER :: linres_section 01470 TYPE(cp_fm_p_type), DIMENSION(:), 01471 POINTER :: vec 01472 INTEGER, INTENT(IN) :: ivec 01473 CHARACTER(LEN=*) :: tag 01474 INTEGER, INTENT(IN), OPTIONAL :: ind 01475 TYPE(cp_error_type), INTENT(INOUT) :: error 01476 01477 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_write_restart', 01478 routineP = moduleN//':'//routineN 01479 01480 CHARACTER(LEN=default_path_length) :: filename 01481 CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status 01482 INTEGER :: handle, i, i_block, ia, ie, 01483 ispin, istat, j, max_block, 01484 nao, nmo, nspins, 01485 output_unit, rst_unit 01486 LOGICAL :: failure 01487 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 01488 TYPE(cp_fm_type), POINTER :: mo_coeff 01489 TYPE(cp_logger_type), POINTER :: logger 01490 TYPE(cp_para_env_type), POINTER :: para_env 01491 TYPE(mo_set_p_type), DIMENSION(:), 01492 POINTER :: mos 01493 TYPE(section_vals_type), POINTER :: print_key 01494 01495 failure = .FALSE. 01496 01497 NULLIFY(logger, mo_coeff, mos, para_env, print_key, vecbuffer) 01498 01499 CALL timeset(routineN,handle) 01500 01501 logger => cp_error_get_logger(error) 01502 01503 IF(BTEST(cp_print_key_should_output(logger%iter_info,linres_section,"PRINT%RESTART",& 01504 used_print_key=print_key,error=error),& 01505 cp_p_file)) THEN 01506 01507 output_unit = cp_print_key_unit_nr(logger,linres_section,& 01508 "PRINT%PROGRAM_RUN_INFO",extension=".Log",error=error) 01509 01510 CALL get_qs_env(qs_env=qs_env, & 01511 mos=mos, & 01512 para_env=para_env, & 01513 error=error) 01514 01515 nspins = SIZE(mos) 01516 01517 my_status = "REPLACE" 01518 my_pos= "REWIND" 01519 CALL XSTRING(tag,ia,ie) 01520 IF(PRESENT(ind)) THEN 01521 my_middle = "RESTART-"//tag(ia:ie)//TRIM(ADJUSTL(cp_to_string(ivec))) 01522 ELSE 01523 my_middle = "RESTART-"//tag(ia:ie) 01524 IF(ivec>1) THEN 01525 my_status="OLD" 01526 my_pos="APPEND" 01527 END IF 01528 END IF 01529 rst_unit = cp_print_key_unit_nr(logger,linres_section,"PRINT%RESTART",& 01530 extension=".lr",middle_name=TRIM(my_middle),file_status=TRIM(my_status),& 01531 file_position=TRIM(my_pos),file_action="WRITE",file_form="UNFORMATTED",& 01532 error=error) 01533 01534 01535 filename = cp_print_key_generate_filename(logger,print_key,& 01536 extension=".lr",middle_name=TRIM(my_middle), my_local=.FALSE.,error=error) 01537 01538 IF(output_unit>0) THEN 01539 WRITE (UNIT=output_unit,FMT="(/,T10,A,A,/)")& 01540 "Linres response functions are written to restart file",& 01541 TRIM(filename) 01542 END IF 01543 01544 ! 01545 ! write data to file 01546 ! use the scalapack block size as a default for buffering columns 01547 CALL get_mo_set(mos(1)%mo_set,mo_coeff=mo_coeff) 01548 CALL cp_fm_get_info(mo_coeff,nrow_global=nao,ncol_block=max_block,error=error) 01549 ALLOCATE(vecbuffer(nao,max_block),STAT=istat) 01550 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01551 01552 IF(PRESENT(ind)) THEN 01553 IF(rst_unit>0) WRITE(rst_unit) ind,ivec,nspins,nao 01554 ELSE 01555 IF(rst_unit>0) WRITE(rst_unit) ivec,nspins,nao 01556 END IF 01557 01558 DO ispin=1,nspins 01559 CALL cp_fm_get_info(vec(ispin)%matrix,ncol_global=nmo,error=error) 01560 01561 IF(rst_unit>0) WRITE(rst_unit) nmo 01562 01563 DO i=1,nmo,MAX(max_block,1) 01564 i_block=MIN(max_block,nmo-i+1) 01565 CALL cp_fm_get_submatrix(vec(ispin)%matrix,vecbuffer,1,i,nao,i_block,error=error) 01566 ! doing this in one write would increase efficiency, but breaks RESTART compatibility. 01567 ! to old ones, and in cases where max_block is different between runs, as might happen during 01568 ! restarts with a different number of CPUs 01569 DO j=1,i_block 01570 IF(rst_unit>0) WRITE (rst_unit) vecbuffer(1:nao,j) 01571 ENDDO 01572 ENDDO 01573 ENDDO 01574 01575 DEALLOCATE (vecbuffer,STAT=istat) 01576 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01577 01578 CALL cp_print_key_finished_output(rst_unit,logger,linres_section,& 01579 "PRINT%RESTART",error=error) 01580 ENDIF 01581 01582 CALL timestop(handle) 01583 01584 END SUBROUTINE linres_write_restart 01585 01586 SUBROUTINE linres_init_write_restart(linres_section,error) 01587 TYPE(section_vals_type), POINTER :: linres_section 01588 TYPE(cp_error_type), INTENT(INOUT) :: error 01589 01590 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_init_write_restart', 01591 routineP = moduleN//':'//routineN 01592 01593 INTEGER :: handle, output_unit, rst_unit 01594 LOGICAL :: failure 01595 TYPE(cp_logger_type), POINTER :: logger 01596 01597 failure = .FALSE. 01598 01599 01600 CALL timeset(routineN,handle) 01601 01602 logger => cp_error_get_logger(error) 01603 01604 IF(BTEST(cp_print_key_should_output(logger%iter_info,linres_section,"PRINT%RESTART",error=error),& 01605 cp_p_file)) THEN 01606 01607 output_unit = cp_print_key_unit_nr(logger,linres_section,& 01608 "PRINT%PROGRAM_RUN_INFO",extension=".Log",error=error) 01609 01610 rst_unit = cp_print_key_unit_nr(logger,linres_section,"PRINT%RESTART",& 01611 extension=".lr",middle_name="RESTART",file_status="REPLACE",& 01612 file_action="WRITE",file_position="APPEND",file_form="UNFORMATTED",& 01613 error=error) 01614 01615 ! may write some infos here about the response... 01616 01617 CALL cp_print_key_finished_output(rst_unit,logger,linres_section,& 01618 "PRINT%RESTART",error=error) 01619 ENDIF 01620 01621 CALL timestop(handle) 01622 01623 END SUBROUTINE linres_init_write_restart 01624 01625 01626 SUBROUTINE linres_read_restart(qs_env,linres_section,vec,ivec,tag,ind,error) 01627 TYPE(qs_environment_type), POINTER :: qs_env 01628 TYPE(section_vals_type), POINTER :: linres_section 01629 TYPE(cp_fm_p_type), DIMENSION(:), 01630 POINTER :: vec 01631 INTEGER, INTENT(IN) :: ivec 01632 CHARACTER(LEN=*) :: tag 01633 INTEGER, INTENT(INOUT), OPTIONAL :: ind 01634 TYPE(cp_error_type), INTENT(INOUT) :: error 01635 01636 CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_read_restart', 01637 routineP = moduleN//':'//routineN 01638 01639 CHARACTER(LEN=default_path_length) :: filename 01640 CHARACTER(LEN=default_string_length) :: my_middle 01641 INTEGER :: group, handle, i, i_block, ia, ie, iostat, ispin, istat, iv, 01642 iv1, ivec_tmp, j, max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, 01643 nspins, nspins_tmp, output_unit, rst_unit, source 01644 LOGICAL :: failure, file_exists 01645 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 01646 TYPE(cp_fm_type), POINTER :: mo_coeff 01647 TYPE(cp_logger_type), POINTER :: logger 01648 TYPE(cp_para_env_type), POINTER :: para_env 01649 TYPE(mo_set_p_type), DIMENSION(:), 01650 POINTER :: mos 01651 TYPE(section_vals_type), POINTER :: print_key 01652 01653 failure = .FALSE. 01654 file_exists = .FALSE. 01655 01656 CALL timeset(routineN,handle) 01657 01658 NULLIFY(mos,para_env,logger,print_key,vecbuffer) 01659 logger => cp_error_get_logger(error) 01660 01661 output_unit = cp_print_key_unit_nr(logger,linres_section,& 01662 "PRINT%PROGRAM_RUN_INFO", extension=".Log",error=error) 01663 01664 CALL get_qs_env(qs_env=qs_env, & 01665 para_env=para_env,& 01666 mos=mos,& 01667 error=error) 01668 01669 nspins = SIZE(mos) 01670 group = para_env%group 01671 source = para_env%source!ionode??? 01672 01673 rst_unit = -1 01674 IF(para_env%ionode) THEN 01675 CALL section_vals_val_get(linres_section,"WFN_RESTART_FILE_NAME",& 01676 n_rep_val=n_rep_val,error=error) 01677 01678 CALL XSTRING(tag,ia,ie) 01679 IF(PRESENT(ind)) THEN 01680 my_middle = "RESTART-"//tag(ia:ie)//TRIM(ADJUSTL(cp_to_string(ivec))) 01681 ELSE 01682 my_middle = "RESTART-"//tag(ia:ie) 01683 END IF 01684 01685 IF(n_rep_val>0) THEN 01686 CALL section_vals_val_get(linres_section,"WFN_RESTART_FILE_NAME",c_val=filename,error=error) 01687 CALL xstring (filename, ia, ie ) 01688 filename = filename(ia:ie)//TRIM(my_middle)//".lr" 01689 ELSE 01690 ! try to read from the filename that is generated automatically from the printkey 01691 print_key => section_vals_get_subs_vals(linres_section,"PRINT%RESTART",error=error) 01692 filename = cp_print_key_generate_filename(logger,print_key, & 01693 extension=".lr",middle_name=TRIM(my_middle),my_local=.FALSE., error=error) 01694 ENDIF 01695 INQUIRE(FILE=filename,exist=file_exists) 01696 ! 01697 ! open file 01698 IF(file_exists) THEN 01699 CALL open_file(file_name=TRIM(filename),& 01700 file_action="READ",& 01701 file_form="UNFORMATTED",& 01702 file_position="REWIND",& 01703 file_status="OLD",& 01704 unit_number=rst_unit) 01705 01706 IF(output_unit > 0) WRITE(output_unit,"(/,T20,A,A)")& 01707 "Read response wavefunction from restart file ", TRIM(filename) 01708 ELSE 01709 IF(output_unit > 0) WRITE(output_unit,"(/,T10,A)")& 01710 "Restart file not available filename=<"//TRIM(filename)//'>' 01711 ENDIF 01712 ENDIF 01713 01714 CALL mp_bcast(file_exists,source,group) 01715 01716 IF(file_exists) THEN 01717 01718 CALL get_mo_set(mos(1)%mo_set,mo_coeff=mo_coeff) 01719 CALL cp_fm_get_info(mo_coeff,nrow_global=nao,ncol_block=max_block,error=error) 01720 01721 ALLOCATE(vecbuffer(nao,max_block),STAT=istat) 01722 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01723 ! 01724 ! read headers 01725 IF(PRESENT(ind)) THEN 01726 iv1 = ivec 01727 ELSE 01728 iv1=1 01729 END IF 01730 DO iv = iv1,ivec 01731 01732 IF(PRESENT(ind)) THEN 01733 IF(rst_unit>0) READ(rst_unit,IOSTAT=iostat) ind, ivec_tmp,nspins_tmp,nao_tmp 01734 CALL mp_bcast(ind,source,group) 01735 ELSE 01736 IF(rst_unit>0) READ(rst_unit,IOSTAT=iostat) ivec_tmp,nspins_tmp,nao_tmp 01737 END IF 01738 CALL mp_bcast(iostat,source,group) 01739 IF(iostat.NE.0) EXIT 01740 CALL mp_bcast(ivec_tmp,source,group) 01741 CALL mp_bcast(nspins_tmp,source,group) 01742 CALL mp_bcast(nao_tmp,source,group) 01743 01744 ! check that the number nao, nmo and nspins are 01745 ! the same as in the current mos 01746 IF(nspins_tmp.NE.nspins) CALL stop_program(routineN,moduleN,__LINE__,& 01747 "nspins not consistent") 01748 IF(nao_tmp .NE.nao ) CALL stop_program(routineN,moduleN,__LINE__,& 01749 "nao not consistent") 01750 ! 01751 DO ispin = 1,nspins 01752 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 01753 CALL cp_fm_get_info(mo_coeff,ncol_global=nmo,error=error) 01754 ! 01755 IF(rst_unit>0) READ(rst_unit) nmo_tmp 01756 CALL mp_bcast(nmo_tmp,source,group) 01757 IF(nmo_tmp.NE.nmo) CALL stop_program(routineN,moduleN,__LINE__,& 01758 "nmo not consistent") 01759 ! 01760 ! read the response 01761 DO i=1,nmo,MAX(max_block,1) 01762 i_block=MIN(max_block,nmo-i+1) 01763 DO j=1,i_block 01764 IF(rst_unit>0) READ(rst_unit) vecbuffer(1:nao,j) 01765 ENDDO 01766 IF(iv.EQ.ivec_tmp) THEN 01767 CALL mp_bcast(vecbuffer,source,group) 01768 CALL cp_fm_set_submatrix(vec(ispin)%matrix,vecbuffer,1,i,nao,i_block,error=error) 01769 ENDIF 01770 ENDDO 01771 ENDDO 01772 IF(ivec.EQ.ivec_tmp) EXIT 01773 ENDDO 01774 01775 IF(iostat.NE.0) THEN 01776 IF(output_unit > 0) WRITE(output_unit,"(/,T10,A,A,/)") & 01777 "Restart file: didnt find ",TRIM(filename) 01778 ENDIF 01779 01780 DEALLOCATE(vecbuffer,STAT=istat) 01781 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01782 01783 ENDIF 01784 01785 IF(para_env%ionode) THEN 01786 IF(file_exists) CALL close_file(unit_number=rst_unit) 01787 ENDIF 01788 01789 CALL timestop(handle) 01790 01791 END SUBROUTINE linres_read_restart 01792 01793 ! ***************************************************************************** 01794 01795 SUBROUTINE check_p_env_init(p_env,linres_control, nspins, error) 01796 ! 01797 TYPE(qs_p_env_type), POINTER :: p_env 01798 TYPE(linres_control_type), POINTER :: linres_control 01799 INTEGER, INTENT(IN) :: nspins 01800 TYPE(cp_error_type), INTENT(INOUT) :: error 01801 01802 CHARACTER(LEN=*), PARAMETER :: routineN = 'check_p_env_init', 01803 routineP = moduleN//':'//routineN 01804 01805 INTEGER :: ispin, ncol, nrow 01806 LOGICAL :: failure 01807 01808 failure = .FALSE. 01809 01810 p_env%iter = 0 01811 p_env%only_energy = .FALSE. 01812 p_env%ls_count=0 01813 01814 p_env%ls_pos = 0.0_dp 01815 p_env%ls_energy = 0.0_dp 01816 p_env%ls_grad = 0.0_dp 01817 p_env%gnorm_old = 1.0_dp 01818 01819 IF(linres_control%preconditioner_type /= ot_precond_none) THEN 01820 CPPrecondition(ASSOCIATED(p_env%preconditioner),cp_failure_level,routineP,error,failure) 01821 DO ispin = 1,nspins 01822 CALL cp_fm_get_info(p_env%PS_psi0(ispin)%matrix,nrow_global=nrow,ncol_global=ncol,error=error) 01823 CPPrecondition(nrow==p_env%n_ao(ispin),cp_failure_level,routineP,error,failure) 01824 CPPrecondition(ncol==p_env%n_mo(ispin),cp_failure_level,routineP,error,failure) 01825 ENDDO 01826 ENDIF 01827 01828 END SUBROUTINE check_p_env_init 01829 01830 END MODULE qs_linres_methods
1.7.3