|
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 ! ***************************************************************************** 00016 MODULE ep_methods 00017 USE atomic_kind_types, ONLY: atomic_kind_type,& 00018 get_atomic_kind_set 00019 USE cp_array_utils, ONLY: cp_2d_r_p_type 00020 USE cp_control_types, ONLY: dft_control_type 00021 USE cp_dbcsr_interface, ONLY: cp_dbcsr_copy,& 00022 cp_dbcsr_init,& 00023 cp_dbcsr_scale,& 00024 cp_dbcsr_set 00025 USE cp_dbcsr_operations, ONLY: cp_dbcsr_allocate_matrix_set,& 00026 cp_dbcsr_deallocate_matrix_set,& 00027 cp_dbcsr_plus_fm_fm_t,& 00028 cp_dbcsr_sm_fm_multiply 00029 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00030 USE cp_files, ONLY: close_file,& 00031 open_file 00032 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm,& 00033 cp_fm_scale_and_add,& 00034 cp_fm_symm,& 00035 cp_fm_trace 00036 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 00037 cp_fm_pool_type,& 00038 fm_pool_create_fm,& 00039 fm_pool_give_back_fm,& 00040 fm_pools_create_fm_vect,& 00041 fm_pools_give_back_fm_vect 00042 USE cp_fm_types, ONLY: cp_fm_get_info,& 00043 cp_fm_get_submatrix,& 00044 cp_fm_p_type,& 00045 cp_fm_set_all,& 00046 cp_fm_set_element,& 00047 cp_fm_set_submatrix,& 00048 cp_fm_type,& 00049 cp_fm_write 00050 USE cp_fm_vect, ONLY: cp_fm_vect_set_all,& 00051 cp_fm_vect_write 00052 USE cp_output_handling, ONLY: & 00053 cp_add_iter_level, cp_iter_string, cp_iterate, cp_p_file, cp_p_store, & 00054 cp_print_key_finished_output, cp_print_key_log, & 00055 cp_print_key_should_output, cp_print_key_unit_nr, cp_rm_iter_level 00056 USE cp_para_env, ONLY: cp_para_env_release,& 00057 cp_para_env_retain 00058 USE cp_para_types, ONLY: cp_blacs_env_type,& 00059 cp_para_env_type 00060 USE cp_subsys_types, ONLY: cp_subsys_get,& 00061 cp_subsys_type 00062 USE ep_qs_methods, ONLY: ep_qs_create 00063 USE ep_qs_types, ONLY: ep_qs_release,& 00064 ep_qs_type 00065 USE ep_types, ONLY: ep_energy_zero,& 00066 ep_env_get,& 00067 ep_env_type,& 00068 ep_envs_get_ep_env,& 00069 ep_force_create,& 00070 ep_force_zero 00071 USE f77_interface, ONLY: calc_energy,& 00072 calc_energy_force,& 00073 create_force_env,& 00074 f_env_add_defaults,& 00075 f_env_rm_defaults,& 00076 f_env_type 00077 USE force_env_types, ONLY: force_env_get 00078 USE global_types, ONLY: global_environment_type 00079 USE input_constants, ONLY: do_ep,& 00080 do_qs,& 00081 ot_precond_full_kinetic,& 00082 ot_precond_solver_default,& 00083 wfi_use_prev_wf_method_nr 00084 USE input_cp2k, ONLY: create_cp2k_input_reading,& 00085 empty_initial_variables 00086 USE input_section_types, ONLY: section_get_ivals,& 00087 section_vals_get_subs_vals,& 00088 section_vals_release,& 00089 section_vals_retain,& 00090 section_vals_type,& 00091 section_vals_val_get,& 00092 section_vals_val_set,& 00093 section_vals_write 00094 USE kahan_sum, ONLY: accurate_sum 00095 USE kinds, ONLY: default_path_length,& 00096 default_string_length,& 00097 dp 00098 USE machine, ONLY: m_flush 00099 USE message_passing, ONLY: mp_bcast,& 00100 mp_max,& 00101 mp_sum 00102 USE particle_list_types, ONLY: particle_list_type 00103 USE particle_types, ONLY: particle_type 00104 USE preconditioner, ONLY: apply_preconditioner,& 00105 make_preconditioner 00106 USE preconditioner_types, ONLY: destroy_preconditioner,& 00107 init_preconditioner 00108 USE qs_charges_types, ONLY: qs_charges_type 00109 USE qs_core_energies, ONLY: calculate_ecore_overlap,& 00110 calculate_ecore_self 00111 USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix 00112 USE qs_energy_types, ONLY: qs_energy_type 00113 USE qs_environment_methods, ONLY: qs_env_update_s_mstruct 00114 USE qs_environment_types, ONLY: get_qs_env,& 00115 qs_env_retain,& 00116 qs_environment_type,& 00117 set_qs_env 00118 USE qs_force, ONLY: write_forces 00119 USE qs_force_types, ONLY: allocate_qs_force,& 00120 qs_force_type,& 00121 zero_qs_force 00122 USE qs_ks_methods, ONLY: qs_ks_create,& 00123 qs_ks_did_change,& 00124 qs_ks_update_qs_env 00125 USE qs_ks_types, ONLY: qs_ks_env_type,& 00126 qs_ks_release 00127 USE qs_matrix_pools, ONLY: mpools_get,& 00128 qs_matrix_pools_type 00129 USE qs_mo_types, ONLY: get_mo_set,& 00130 init_mo_set,& 00131 mo_set_p_type 00132 USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists 00133 USE qs_p_env_methods, ONLY: p_env_create,& 00134 p_env_did_change,& 00135 p_env_psi0_changed,& 00136 p_op_l1,& 00137 p_op_l2_fawzi,& 00138 p_postortho,& 00139 p_preortho 00140 USE qs_p_env_types, ONLY: qs_p_env_type 00141 USE qs_p_sparse_psi, ONLY: p_proj_create,& 00142 qs_p_projection_p_type,& 00143 qs_p_projection_type 00144 USE qs_rho_methods, ONLY: qs_rho_update_rho 00145 USE qs_rho_types, ONLY: qs_rho_type 00146 USE qs_wf_history_methods, ONLY: wfi_change_memory_depth,& 00147 wfi_create,& 00148 wfs_duplicate_snapshot 00149 USE qs_wf_history_types, ONLY: qs_wf_history_type,& 00150 qs_wf_snapshot_type,& 00151 wfi_get_snapshot,& 00152 wfi_release,& 00153 wfs_release 00154 USE realspace_grid_cube, ONLY: pw_to_cube 00155 USE replica_types, ONLY: rep_env_calc_e_f,& 00156 rep_env_create,& 00157 rep_env_local_index,& 00158 replica_env_type 00159 USE scf_control_types, ONLY: scf_control_type 00160 USE timings, ONLY: timeset,& 00161 timestop 00162 #include "cp_common_uses.h" 00163 00164 IMPLICIT NONE 00165 PRIVATE 00166 00167 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00168 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ep_methods' 00169 INTEGER, PRIVATE, SAVE :: last_ep_env_id=0 00170 00171 PUBLIC :: ep_env_init_low, ep_env_transfer_psi0, & 00172 ep_env_localize_matrix, ep_env_calc_e_f_low 00173 00174 CONTAINS 00175 00176 ! ***************************************************************************** 00182 SUBROUTINE print_qs_energies(qs_env,output_unit,error) 00183 TYPE(qs_environment_type), POINTER :: qs_env 00184 INTEGER, INTENT(IN) :: output_unit 00185 TYPE(cp_error_type), INTENT(inout) :: error 00186 00187 CHARACTER(len=*), PARAMETER :: routineN = 'print_qs_energies', 00188 routineP = moduleN//':'//routineN 00189 00190 INTEGER :: handle 00191 LOGICAL :: failure 00192 TYPE(dft_control_type), POINTER :: dft_control 00193 TYPE(qs_charges_type), POINTER :: qs_charges 00194 TYPE(qs_energy_type), POINTER :: energy 00195 TYPE(qs_rho_type), POINTER :: rho 00196 00197 failure=.FALSE. 00198 00199 CALL timeset(routineN,handle) 00200 IF (output_unit>0) THEN 00201 NULLIFY(rho,energy,qs_charges,dft_control) 00202 CALL get_qs_env(qs_env,rho=rho,energy=energy,dft_control=dft_control,& 00203 qs_charges=qs_charges,error=error) 00204 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T61,F20.10))")& 00205 "Total electronic density (r-space): ",& 00206 accurate_sum(rho%tot_rho_r),& 00207 "Total core charge density (r-space):",& 00208 qs_charges%total_rho_core_rspace 00209 WRITE (UNIT=output_unit,FMT="(T3,A,T61,F20.10)")& 00210 "Total charge density (r-space): ",& 00211 accurate_sum(rho%tot_rho_r)+& 00212 qs_charges%total_rho_core_rspace,& 00213 "Total charge density (g-space): ",qs_charges%total_rho_gspace 00214 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 00215 "Overlap energy of the core charge distribution:",energy%core_overlap,& 00216 "Self energy of the core charge distribution: ",energy%core_self,& 00217 "Core Hamiltonian energy: ",energy%core,& 00218 "Hartree energy: ",energy%hartree,& 00219 "Exchange-correlation energy: ",energy%exc 00220 IF (energy%e_hartree /= 0.0_dp) & 00221 WRITE (UNIT=output_unit,FMT="(T3,A,/,T3,A,T56,F25.14)")& 00222 "Coulomb Electron-Electron Interaction Energy ",& 00223 "- Already included in the total Hartree term ",energy%e_hartree 00224 IF (dft_control%qs_control%mulliken_restraint) THEN 00225 WRITE (UNIT=output_unit,FMT="(T3,A,T61,F20.10)")& 00226 "Mulliken restraint energy: ",energy%mulliken 00227 END IF 00228 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 00229 "Total energy: ",energy%total 00230 CALL m_flush(output_unit) 00231 END IF 00232 CALL timestop(handle) 00233 END SUBROUTINE print_qs_energies 00234 00235 ! ***************************************************************************** 00244 SUBROUTINE calc_c0_tilde(ep_env,mol_pos,new_c,imol,error) 00245 TYPE(ep_env_type), POINTER :: ep_env 00246 REAL(kind=dp), DIMENSION(:), INTENT(in) :: mol_pos 00247 TYPE(cp_2d_r_p_type), DIMENSION(:), 00248 INTENT(inout) :: new_c 00249 INTEGER, INTENT(in) :: imol 00250 TYPE(cp_error_type), INTENT(inout) :: error 00251 00252 INTEGER :: icol, irow, ispin 00253 REAL(dp), DIMENSION(:), POINTER :: coeff_standard1d 00254 00255 ! legge dall'input 00256 00257 CALL section_vals_val_get(ep_env%root_section,"FORCE_EVAL%EP%START_COEFFS",& 00258 r_vals=coeff_standard1d,error=error) 00259 00260 ! trova rotazione usando la posizione mol_pos (coordinate xyz xyz xyz) 00261 ! che porta la molecola nella posizione standard 00262 00263 ! ruota gli orbitali standard per avere gli orbitali della molecola 00264 ! attuale in new_c 00265 DO ispin=1,ep_env%nspins 00266 DO icol=1,ep_env%sub_nmo(ispin) 00267 DO irow=1,ep_env%sub_nao(ispin) 00268 new_c(ispin)%array(icol,irow)=0.0_dp 00269 END DO 00270 END DO 00271 END DO 00272 END SUBROUTINE calc_c0_tilde 00273 00274 ! ***************************************************************************** 00288 SUBROUTINE ep_env_init_low(ep_env_id,ierr) 00289 INTEGER, INTENT(in) :: ep_env_id 00290 INTEGER, INTENT(out) :: ierr 00291 00292 CHARACTER(len=*), PARAMETER :: routineN = 'ep_env_init_low', 00293 routineP = moduleN//':'//routineN 00294 00295 CHARACTER(len=default_path_length) :: comp_path 00296 CHARACTER(LEN=default_string_length) :: project_name 00297 INTEGER :: at_per_mol, handle, i, iat, 00298 imol, ispin, nat, nmol, stat, 00299 unit_nr 00300 INTEGER, ALLOCATABLE, DIMENSION(:) :: atoms 00301 LOGICAL :: failure 00302 TYPE(cp_blacs_env_type), POINTER :: blacs_env 00303 TYPE(cp_error_type) :: error, new_error 00304 TYPE(cp_fm_pool_p_type), DIMENSION(:), 00305 POINTER :: ao_mo_fm_pools 00306 TYPE(cp_logger_type), POINTER :: logger 00307 TYPE(cp_para_env_type), POINTER :: para_env 00308 TYPE(dft_control_type), POINTER :: dft_control 00309 TYPE(ep_env_type), POINTER :: ep_env 00310 TYPE(ep_qs_type), POINTER :: ep_qs_env 00311 TYPE(f_env_type), POINTER :: f_env 00312 TYPE(global_environment_type), POINTER :: globenv 00313 TYPE(mo_set_p_type), DIMENSION(:), 00314 POINTER :: full_mos, sub_mos 00315 TYPE(particle_type), DIMENSION(:), 00316 POINTER :: particle_set 00317 TYPE(qs_environment_type), POINTER :: sub_qs_env 00318 TYPE(qs_matrix_pools_type), POINTER :: mpools 00319 TYPE(qs_p_projection_p_type), 00320 DIMENSION(:), POINTER :: sub_proj 00321 TYPE(section_vals_type), POINTER :: comp_input, input 00322 00323 ! TYPE(qs_p_projection_p_type), & 00324 ! DIMENSION(:), INTENT(in) :: projections 00325 ! TYPE(cp_para_env_type), POINTER :: para_env 00326 00327 CALL timeset(routineN,handle) 00328 failure=.FALSE. 00329 ep_env => ep_envs_get_ep_env(ep_env_id) 00330 NULLIFY(particle_set) 00331 CALL cp_error_init(error) 00332 CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 00333 IF (.NOT.failure) THEN 00334 input => ep_env%root_section 00335 para_env => ep_env%para_env 00336 logger => cp_error_get_logger(error) 00337 CALL cp_log(logger, level=cp_note_level, fromWhere=routineP ,& 00338 message=routineP//" creating main_qs_env", local=.FALSE.) 00339 00340 CALL cp_print_key_log(logger,input,"FORCE_EVAL%EP%PRINT%RUN_INFO",& 00341 extension=".epLog",message=routineP//" creating main_qs_env",error=error) 00342 ! create_main_qs_env 00343 CALL section_vals_val_set(input,"FORCE_EVAL%METHOD",& 00344 i_val=do_qs,error=error) 00345 CALL section_vals_val_get(input,"GLOBAL%PROJECT_NAME",& 00346 c_val=project_name,error=error) 00347 CALL section_vals_val_set(input,"GLOBAL%PROJECT_NAME",& 00348 c_val=TRIM(project_name)//"_mainQS",error=error) 00349 CALL open_file(file_name=TRIM(project_name)//"_mainQS.inp",& 00350 file_status="UNKNOWN",file_form="FORMATTED",file_action="WRITE",& 00351 unit_number=unit_nr) 00352 CALL section_vals_write(input,unit_nr=unit_nr,hide_root=.TRUE.,& 00353 error=error) 00354 CALL close_file(unit_nr) 00355 CALL section_vals_val_set(input,"GLOBAL%PROJECT_NAME",& 00356 c_val=project_name,error=error) 00357 CALL section_vals_val_set(input,"FORCE_EVAL%METHOD",& 00358 i_val=do_ep,error=error) 00359 00360 CALL create_force_env(ep_env%f_env_id,TRIM(project_name)//"_mainQS.inp",& 00361 TRIM(project_name)//"_mainQS.out",para_env%group,ierr=ierr) 00362 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 00363 globenv => ep_env%globenv 00364 NULLIFY(f_env) 00365 CALL f_env_add_defaults(f_env_id=ep_env%f_env_id,f_env=f_env,& 00366 new_error=new_error, failure=failure) 00367 CALL force_env_get(f_env%force_env,qs_env=ep_env%main_qs_env,& 00368 globenv=globenv,error=new_error) 00369 CALL qs_env_retain(ep_env%main_qs_env,error=error) 00370 CALL cp_print_key_log(logger,input,"FORCE_EVAL%EP%PRINT%RUN_INFO",& 00371 extension=".epLog",message=routineP//" created main_qs_env "//& 00372 TRIM(ADJUSTR(cp_to_string(ep_env%main_qs_env%id_nr))),error=error) 00373 00374 !- ! reduce memory usage, get rid of an input structure 00375 !- CALL section_vals_retain(root_section,error=error) 00376 !- CALL section_vals_release(ep_env%root_section,error=error) 00377 !- ep_env%root_section => root_section 00378 !- input => ep_env%root_section 00379 CALL section_vals_release(ep_env%input,error=error) 00380 ep_env%input => section_vals_get_subs_vals(ep_env%root_section,& 00381 "FORCE_EVAL%EP",error=error) 00382 CALL section_vals_retain(ep_env%input,error=error) 00383 !- CALL globenv_retain(globenv,error=error) 00384 !- CALL globenv_release(ep_env%globenv,error=error) 00385 !- ep_env%globenv => globenv 00386 00387 ! *** initial initializations (atom position dependent) 00388 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00389 extension=".epLog",message=routineP//" initial setup of main_qs_env",& 00390 error=error) 00391 00392 CALL ep_env_finish_qs_init(ep_env%main_qs_env, globenv, error) 00393 00394 CALL get_qs_env(ep_env%main_qs_env,dft_control=dft_control,& 00395 blacs_env=blacs_env,para_env=para_env,error=error) 00396 00397 CALL cp_para_env_retain(para_env,error=error) 00398 CALL cp_para_env_release(ep_env%para_env,error=error) 00399 ep_env%para_env => para_env 00400 00401 ep_env%nspins=dft_control%nspins 00402 ALLOCATE(ep_env%full_nao(ep_env%nspins),ep_env%full_nmo(ep_env%nspins),& 00403 stat=stat) 00404 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00405 NULLIFY(full_mos) 00406 CALL get_qs_env(ep_env%main_qs_env,mos=full_mos,& 00407 particle_set=particle_set,mpools=mpools,error=new_error) 00408 DO ispin=1,ep_env%nspins 00409 CALL get_mo_set(full_mos(ispin)%mo_set,nao=ep_env%full_nao(ispin),& 00410 nmo=ep_env%full_nmo(ispin)) 00411 END DO 00412 00413 ! add ep_qs_env 00414 NULLIFY(ep_qs_env) 00415 CALL ep_qs_create(ep_qs_env,ep_env%main_qs_env,error=error) 00416 ! CALL set_qs_env(ep_env%main_qs_env,ep_qs_env=ep_qs_env,error=error) 00417 CALL ep_qs_release(ep_qs_env,error=error) 00418 ! 00419 CALL f_env_rm_defaults(f_env,new_error,ierr) 00420 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 00421 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00422 extension=".epLog",message=routineP//& 00423 " initial setup of main_qs_env finished",error=error) 00424 00425 ! create mol_envs 00426 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00427 extension=".epLog",message=routineP//& 00428 " creating mol_envs",error=error) 00429 CALL cp_get_natom(ep_env%f_env_id,nat,ierr) 00430 CALL section_vals_val_get(input,"FORCE_EVAL%EP%AT_PER_MOL",& 00431 i_val=at_per_mol,error=error) 00432 nmol=nat/at_per_mol 00433 CPPrecondition(MODULO(nat,at_per_mol)==0,cp_failure_level,routineP,error,failure) 00434 ep_env%nmol=nmol 00435 ep_env%nat=nat 00436 ep_env%nat_per_mol=at_per_mol 00437 CALL ep_force_create(ep_env%force,nat=nat,error=error) 00438 CALL section_vals_val_get(input,"FORCE_EVAL%EP%COMP_INPUT",& 00439 c_val=comp_path,error=error) 00440 comp_input => create_cp2k_input_reading(comp_path, initial_variables=empty_initial_variables, & 00441 para_env=para_env,error=error) 00442 CALL rep_env_create(ep_env%mol_envs, para_env=para_env, input=comp_input,& 00443 nrep=nmol, prep=1,keep_wf_history=.TRUE.,error=error) 00444 CALL cp_assert(ASSOCIATED(ep_env%mol_envs),cp_failure_level,& 00445 cp_assertion_failed,routineP,"ep not implemented for a number of"//& 00446 " processors bigger than the number of molecules",error,failure) 00447 CALL section_vals_release(comp_input,error=error) 00448 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00449 extension=".epLog",message=routineP//& 00450 " mol_envs created",error=error) 00451 00452 !sub_nao,sub_nmo 00453 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00454 extension=".epLog",message=routineP//& 00455 " mol_envs initial setup",error=error) 00456 ALLOCATE(ep_env%sub_nao(ep_env%nspins),ep_env%sub_nmo(ep_env%nspins),& 00457 stat=stat) 00458 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00459 ep_env%sub_nao=0 00460 ep_env%sub_nmo=0 00461 IF (ASSOCIATED(ep_env%mol_envs)) THEN 00462 NULLIFY(sub_mos,sub_qs_env) 00463 CALL f_env_add_defaults(f_env_id=ep_env%mol_envs%f_env_id,f_env=f_env,& 00464 new_error=new_error, failure=failure) 00465 CALL force_env_get(f_env%force_env,qs_env=sub_qs_env,error=new_error) 00466 CALL get_qs_env(sub_qs_env,mos=sub_mos,error=new_error) 00467 DO ispin=1,ep_env%nspins 00468 CALL get_mo_set(sub_mos(ispin)%mo_set,nao=ep_env%sub_nao(ispin),& 00469 nmo=ep_env%sub_nmo(ispin)) 00470 END DO 00471 00472 CALL f_env_rm_defaults(f_env,new_error,ierr) 00473 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 00474 END IF 00475 CALL mp_max(ep_env%sub_nmo,para_env%group) 00476 CALL mp_max(ep_env%sub_nao,para_env%group) 00477 ! dummy force calc to ensure that everything is allocated for force 00478 ! evaluation 00479 CALL calc_energy_force(ep_env%mol_envs%f_env_id,.TRUE.,ierr) 00480 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 00481 ! guarantee wf & overlap in wf_history 00482 DO imol=1,SIZE(ep_env%mol_envs%local_rep_indices) 00483 ep_env%mol_envs%wf_history(imol)%wf_history%store_wf=.TRUE. 00484 ep_env%mol_envs%wf_history(imol)%wf_history%store_overlap=.TRUE. 00485 CALL wfi_change_memory_depth(& 00486 ep_env%mol_envs%wf_history(imol)%wf_history,& 00487 MAX(1,ep_env%mol_envs%wf_history(imol)%wf_history%memory_depth),& 00488 error=error) 00489 END DO 00490 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00491 extension=".epLog",message=routineP//& 00492 " mol_envs setup fininshed",error=error) 00493 ! create the perturbation environment 00494 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00495 extension=".epLog",message=routineP//& 00496 " creating sub_p_env",error=error) 00497 CALL p_env_create(ep_env%sub_p_env, qs_env=sub_qs_env,& 00498 orthogonal_orbitals=.TRUE.,error=error) 00499 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00500 extension=".epLog",message=routineP//& 00501 " sub_p_env created",error=error) 00502 00503 ! create sub_proj and at2sub 00504 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00505 extension=".epLog",message=routineP//& 00506 " creating projections",error=error) 00507 ALLOCATE(sub_proj(nmol),stat=stat) 00508 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00509 ALLOCATE(ep_env%at2sub(nat),stat=stat) 00510 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00511 ALLOCATE(atoms(at_per_mol),stat=stat) 00512 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00513 iat=0 00514 DO imol=1,nmol 00515 DO i=1,at_per_mol 00516 iat=iat+1 00517 atoms(i)=iat 00518 ep_env%at2sub(iat)=imol 00519 END DO 00520 CALL p_proj_create(sub_proj(imol)%projection, atoms=atoms,& 00521 particle_set=particle_set, error=error) 00522 END DO 00523 DEALLOCATE(atoms,stat=stat) 00524 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 00525 ep_env%sub_proj => sub_proj 00526 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00527 extension=".epLog",message=routineP//& 00528 " created projections",error=error) 00529 00530 ! create the perturbation environment 00531 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00532 extension=".epLog",message=routineP//& 00533 " creating main_p_env",error=error) 00534 CALL p_env_create(ep_env%main_p_env, qs_env=ep_env%main_qs_env,& 00535 orthogonal_orbitals=.FALSE.,error=error) 00536 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00537 extension=".epLog",message=routineP//& 00538 " created main_p_env",error=error) 00539 00540 ! alloc m_pi_Hrho_psi0d 00541 CALL get_qs_env(ep_env%main_qs_env, mpools=mpools,& 00542 error=error) 00543 CALL mpools_get(mpools,ao_mo_fm_pools=ao_mo_fm_pools, error=error) 00544 00545 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,ep_env%m_pi_Hrho_psi0d,& 00546 name="ep_env"//TRIM(ADJUSTL(cp_to_string(ep_env%id_nr)))//& 00547 "%m_pi_Hrho_psi0d",error=error) 00548 ! alloc psi1 00549 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,ep_env%psi1,& 00550 name="ep_env"//TRIM(ADJUSTL(cp_to_string(ep_env%id_nr)))//"%psi1",& 00551 error=error) 00552 00553 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 00554 extension=".epLog",message=routineP//& 00555 " ep initial setup finished",error=error) 00556 END IF 00557 CALL cp_error_check(error,failure) 00558 IF (failure) THEN 00559 ierr=1 00560 ELSE 00561 ierr=0 00562 END IF 00563 CALL cp_error_dealloc_ref(error) 00564 CALL timestop(handle) 00565 END SUBROUTINE ep_env_init_low 00566 00567 ! ***************************************************************************** 00577 SUBROUTINE ep_env_finish_qs_init(qs_env,globenv,error) 00578 TYPE(qs_environment_type), POINTER :: qs_env 00579 TYPE(global_environment_type), POINTER :: globenv 00580 TYPE(cp_error_type), INTENT(inout) :: error 00581 00582 CHARACTER(len=*), PARAMETER :: routineN = 'ep_env_finish_qs_init', 00583 routineP = moduleN//':'//routineN 00584 00585 INTEGER :: handle, ispin 00586 LOGICAL :: failure 00587 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00588 POINTER :: matrix_ks, matrix_s 00589 TYPE(cp_fm_type), POINTER :: mo_coeff 00590 TYPE(cp_para_env_type), POINTER :: para_env 00591 TYPE(dft_control_type), POINTER :: dft_control 00592 TYPE(mo_set_p_type), DIMENSION(:), 00593 POINTER :: mos 00594 TYPE(qs_ks_env_type), POINTER :: ks_env 00595 00596 failure=.FALSE. 00597 NULLIFY(ks_env, mos, mo_coeff, dft_control,matrix_ks,matrix_s) 00598 00599 CALL timeset(routineN,handle) 00600 para_env=>qs_env%para_env 00601 00602 CALL build_qs_neighbor_lists(qs_env,para_env,force_env_section=qs_env%input,error=error) 00603 00604 ! *** Calculate the overlap and the core Hamiltonian integral matrix *** 00605 00606 CALL build_core_hamiltonian_matrix(qs_env=qs_env,calculate_forces=.FALSE.,error=error) 00607 CALL qs_env_update_s_mstruct(qs_env, error=error) 00608 00609 ! *** init the ks_env 00610 00611 CALL get_qs_env(qs_env, ks_env=ks_env,error=error) 00612 IF (.not.ASSOCIATED(ks_env)) THEN 00613 CALL qs_ks_create(ks_env,qs_env=qs_env,error=error) 00614 CALL set_qs_env(qs_env, ks_env=ks_env,error=error) 00615 CALL qs_ks_release(ks_env,error=error) 00616 END IF 00617 00618 ! *** Initializes the MOs *** 00619 CALL get_qs_env(qs_env,mos=mos,error=error) 00620 CPPrecondition(ASSOCIATED(mos),cp_failure_level,routineP,error,failure) 00621 IF (.NOT.failure) THEN 00622 DO ispin=1,SIZE(mos) 00623 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 00624 IF (.NOT.ASSOCIATED(mo_coeff)) THEN 00625 CALL init_mo_set(mos(ispin)%mo_set, & 00626 qs_env%mpools%ao_mo_fm_pools(ispin)%pool,& 00627 name="qs_env"//TRIM(ADJUSTL(cp_to_string(qs_env%id_nr)))//& 00628 "%mo"//TRIM(ADJUSTL(cp_to_string(ispin))),& 00629 error=error) 00630 END IF 00631 END DO 00632 END IF 00633 00634 ! *** Allocate k and put it in the QS environment *** 00635 CALL get_qs_env(qs_env,dft_control=dft_control,matrix_s=matrix_s,error=error) 00636 00637 CALL cp_dbcsr_allocate_matrix_set(matrix_ks,dft_control%nspins,error=error) 00638 DO ispin=1,dft_control%nspins 00639 ALLOCATE (matrix_ks(ispin)%matrix) 00640 CALL cp_dbcsr_init(matrix_ks(ispin)%matrix,error=error) 00641 CALL cp_dbcsr_copy(matrix_ks(ispin)%matrix,matrix_s(1)%matrix,& 00642 "KOHN-SHAM MATRIX-"//ADJUSTL(cp_to_string(ispin)),error=error) 00643 END DO 00644 00645 CALL set_qs_env(qs_env=qs_env,matrix_ks=matrix_ks,error=error) 00646 00647 NULLIFY(qs_env%atprop) 00648 CALL calculate_ecore_self(qs_env,error=error) 00649 CALL calculate_ecore_overlap(qs_env,para_env,.FALSE.,error=error) 00650 00651 CALL timestop(handle) 00652 END SUBROUTINE ep_env_finish_qs_init 00653 00654 ! ***************************************************************************** 00664 SUBROUTINE ep_env_transfer_psi0(ep_env, error) 00665 TYPE(ep_env_type), POINTER :: ep_env 00666 TYPE(cp_error_type), INTENT(inout) :: error 00667 00668 CHARACTER(len=*), PARAMETER :: routineN = 'ep_env_transfer_psi0', 00669 routineP = moduleN//':'//routineN 00670 INTEGER, PARAMETER :: max_blocksize = 100 00671 00672 INTEGER :: blocksize, blocksize2, full_nmo, handle, i, icol, ierr, ispin, 00673 isub, local_idx, my_start_col_full, my_start_col_min, ncol_min, nmo, 00674 nrow_min, nspins, stat, sub_nmo, unit_nr 00675 LOGICAL :: failure, print_matrixes 00676 REAL(kind=dp), ALLOCATABLE, 00677 DIMENSION(:, :) :: sub_mo_matrix, tmp_full 00678 TYPE(cp_error_type) :: new_error, suberror 00679 TYPE(cp_fm_type), POINTER :: orbitals, sub_orbitals 00680 TYPE(cp_logger_type), POINTER :: logger, new_logger 00681 TYPE(f_env_type), POINTER :: f_env 00682 TYPE(mo_set_p_type), DIMENSION(:), 00683 POINTER :: mos 00684 TYPE(qs_environment_type), POINTER :: main_qs_env 00685 TYPE(qs_p_projection_p_type), 00686 DIMENSION(:), POINTER :: sub_proj 00687 TYPE(qs_p_projection_type), POINTER :: proj 00688 TYPE(qs_wf_snapshot_type), POINTER :: snapshot 00689 TYPE(replica_env_type), POINTER :: rep_env 00690 00691 CALL timeset(routineN,handle) 00692 00693 failure=.FALSE. 00694 NULLIFY(sub_orbitals, main_qs_env, sub_proj, logger, orbitals) 00695 logger => cp_error_get_logger(error) 00696 00697 CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 00698 CALL ep_env_get(ep_env,& 00699 sub_proj=sub_proj,main_qs_env=main_qs_env,error=error) 00700 CALL cp_error_check(error,failure) 00701 CPPrecondition(ASSOCIATED(sub_proj),cp_failure_level,routineP,error,failure) 00702 CPPrecondition(ASSOCIATED(main_qs_env),cp_failure_level,routineP,error,failure) 00703 IF (.NOT. failure) THEN 00704 nspins=main_qs_env%dft_control%nspins 00705 CALL get_qs_env(main_qs_env,mos=mos,error=error) 00706 rep_env => ep_env%mol_envs 00707 00708 DO ispin=1,nspins ! transfer also occupation_numbers ? 00709 my_start_col_full=1 00710 sub_nmo=ep_env%sub_nmo(ispin) 00711 full_nmo=ep_env%full_nmo(ispin) 00712 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=orbitals,nmo=nmo) 00713 CALL cp_fm_set_all(orbitals,0.0_dp,error=error) 00714 blocksize=MIN(max_blocksize,sub_nmo) 00715 ALLOCATE(sub_mo_matrix(blocksize,ep_env%sub_nao(ispin)),stat=stat) 00716 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00717 ALLOCATE(tmp_full(blocksize,ep_env%full_nao(ispin)),stat=stat) 00718 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 00719 print_matrixes=BTEST(cp_print_key_should_output(logger%iter_info,& 00720 ep_env%input,"PRINT%EP_MATRIXES/PSI0_BLOCKS",error=error),cp_p_file) 00721 DO isub=1,ep_env%nmol 00722 ! does not work if there are "empty" processors in replica env 00723 CPAssert(ASSOCIATED(rep_env),cp_failure_level,routineP,error,failure) 00724 00725 my_start_col_min=1 00726 proj => ep_env%sub_proj(isub)%projection 00727 blocksize2=blocksize 00728 DO icol=0,sub_nmo-1,blocksize 00729 IF (icol+blocksize>sub_nmo) blocksize2=sub_nmo-icol 00730 IF (rep_env%rep_is_local(isub)) THEN 00731 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id,f_env=f_env,& 00732 new_error=new_error, failure=failure) 00733 local_idx=rep_env_local_index(rep_env,isub,error=error) 00734 CPAssert(local_idx>0,cp_failure_level,routineP,error,failure) 00735 new_logger => cp_error_get_logger(new_error) 00736 snapshot => wfi_get_snapshot(& 00737 rep_env%wf_history(local_idx)%wf_history,& 00738 1,error=new_error) 00739 CPAssert(ASSOCIATED(snapshot),cp_failure_level,routineP,new_error,failure) 00740 CPAssert(ASSOCIATED(snapshot%wf),cp_failure_level,routineP,new_error,failure) 00741 sub_orbitals => snapshot%wf(ispin)%matrix 00742 new_logger%iter_info%iteration(1)=isub 00743 IF (print_matrixes) THEN 00744 CALL cp_error_init(suberror,template_error=new_error) 00745 unit_nr=cp_print_key_unit_nr(new_logger,ep_env%input,& 00746 "PRINT%EP_MATRIXES/PSI0_BLOCKS",extension=".psi0_block",& 00747 ignore_should_output=.TRUE.,error=new_error) 00748 CALL cp_fm_write(sub_orbitals, unit_nr,& 00749 long_description=.TRUE., error=suberror) 00750 CALL cp_error_reset(suberror) 00751 CALL cp_print_key_finished_output(unit_nr,new_logger,ep_env%input,& 00752 "PRINT%EP_MATRIXES/PSI0_BLOCKS",& 00753 ignore_should_output=.TRUE.,error=new_error) 00754 CALL cp_error_dealloc_ref(suberror) 00755 END IF 00756 CALL cp_fm_get_info(sub_orbitals, nrow_global=nrow_min, ncol_global=ncol_min,& 00757 error=error) 00758 CPPostcondition(nrow_min==ep_env%sub_nao(ispin),cp_failure_level,routineP,error,failure) 00759 CPPostcondition(ncol_min==sub_nmo,cp_failure_level,routineP,error,failure) 00760 CALL cp_assert(blocksize2+icol<=sub_nmo+1,& 00761 cp_failure_level,cp_assertion_failed,routineP,& 00762 "icol+blocksize2>sub_nmo"//& 00763 CPSourceFileRef,& 00764 error,failure) 00765 CALL cp_fm_get_submatrix(sub_orbitals,sub_mo_matrix,transpose=.TRUE.,& 00766 start_col=icol+1, n_cols=blocksize2,error=new_error) 00767 CALL f_env_rm_defaults(f_env,new_error,ierr) 00768 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 00769 END IF 00770 CALL mp_bcast(sub_mo_matrix,rep_env%inter_rep_rank(rep_env%replica_owner(isub)),& 00771 rep_env%para_env_inter_rep%group) 00772 00773 ! be smarter and use only a part of the full width? 00774 ! i.e.size(tmp_full,2) could easily reduced to (min(proj_indexes):max(proj_indexes)) 00775 CALL dcopy(SIZE(tmp_full,1)*SIZE(tmp_full,2),0.0_dp,0,tmp_full(1,1),1) 00776 DO i=1,SIZE(proj%proj_indexes) 00777 tmp_full(1:blocksize2,proj%proj_indexes(i))=sub_mo_matrix(1:blocksize2,i) 00778 END DO 00779 00780 CALL cp_assert(icol+my_start_col_full+blocksize2<=ep_env%full_nmo(ispin)+1,& 00781 cp_failure_level,cp_assertion_failed,routineP,& 00782 'icol+my_start_col_full+blocksize2>ep_env%full_nmo(ispin)+1 '//& 00783 CPSourceFileRef,& 00784 error,failure) 00785 CALL cp_fm_set_submatrix(orbitals,tmp_full,& 00786 start_col=icol+my_start_col_full, n_cols=blocksize2, & 00787 transpose=.TRUE.,error=error) 00788 END DO 00789 my_start_col_full=my_start_col_full+sub_nmo 00790 CALL cp_assert(my_start_col_full<=full_nmo+1,& 00791 cp_failure_level,cp_assertion_failed,routineP,& 00792 "my_start_col_full>full_nmo+1, "//& 00793 CPSourceFileRef,& 00794 error,failure) 00795 00796 END DO 00797 CALL cp_assert(my_start_col_full==ep_env%full_nmo(ispin)+1,& 00798 cp_failure_level,cp_assertion_failed,routineP,& 00799 "my_start_col_full/=ep_env%full_nmo(ispin)+1, "//& 00800 CPSourceFileRef,& 00801 error,failure) 00802 00803 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00804 ep_env%input,"PRINT%EP_MATRIXES/PSI0",error=error),cp_p_file)) THEN 00805 CALL cp_error_init(suberror,template_error=error) 00806 unit_nr=cp_print_key_unit_nr(logger,ep_env%input,& 00807 "PRINT%EP_MATRIXES/PSI0",extension=".psi0",& 00808 error=error) 00809 CALL cp_fm_write(orbitals, unit_nr,& 00810 long_description=.TRUE., error=suberror) 00811 CALL cp_error_reset(suberror) 00812 CALL cp_print_key_finished_output(unit_nr,logger,ep_env%input,& 00813 "PRINT%EP_MATRIXES/PSI0",& 00814 error=error) 00815 CALL cp_error_dealloc_ref(suberror) 00816 END IF 00817 END DO 00818 00819 END IF 00820 CALL timestop(handle) 00821 END SUBROUTINE ep_env_transfer_psi0 00822 00823 ! ***************************************************************************** 00839 SUBROUTINE stupid_solve(ep_env,eps_r,error) 00840 TYPE(ep_env_type), POINTER :: ep_env 00841 REAL(KIND=dp), INTENT(in) :: eps_r 00842 TYPE(cp_error_type), INTENT(inout) :: error 00843 00844 CHARACTER(len=*), PARAMETER :: routineN = 'stupid_solve', 00845 routineP = moduleN//':'//routineN 00846 00847 INTEGER :: i, ispin, j, nspins, 00848 output_unit 00849 LOGICAL :: failure 00850 REAL(KIND=dp) :: alpha, Ap_p, Ap_p_s, beta, p_ortho_norm, p_ortho_norm_s, 00851 r_norm2, r_norm2_s, r_z, r_z_new, r_z_s 00852 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00853 POINTER :: matrix_s 00854 TYPE(cp_fm_p_type), DIMENSION(:), 00855 POINTER :: Ap, p, p_ortho, r, z 00856 TYPE(cp_fm_pool_p_type), DIMENSION(:), 00857 POINTER :: ao_mo_fm_pools 00858 TYPE(cp_logger_type), POINTER :: logger 00859 TYPE(dft_control_type), POINTER :: dft_control 00860 TYPE(qs_matrix_pools_type), POINTER :: mpools 00861 00862 failure=.FALSE. 00863 logger => cp_error_get_logger(error) 00864 NULLIFY(p,p_ortho,Ap,r,z,ao_mo_fm_pools,dft_control,matrix_s,mpools) 00865 00866 CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 00867 CPPrecondition(ep_env%ref_count>0,cp_failure_level,routineP,error,failure) 00868 IF (.NOT. failure) THEN 00869 CALL get_qs_env(ep_env%main_qs_env, & 00870 dft_control=dft_control,matrix_s=matrix_s,mpools=mpools,error=error) 00871 00872 CALL mpools_get(mpools,ao_mo_fm_pools=ao_mo_fm_pools, error=error) 00873 nspins=dft_control%nspins 00874 00875 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,p,name="ep_CG_p",& 00876 error=error) 00877 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,p_ortho,& 00878 name="ep_CG_p_ortho",& 00879 error=error) 00880 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,Ap,name="ep_CG_Ap",& 00881 error=error) 00882 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,r,name="ep_CG_r",& 00883 error=error) 00884 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,z,name="ep_CG_z",& 00885 error=error) 00886 CALL cp_error_check(error,failure) 00887 END IF 00888 IF (.not.failure) THEN 00889 CALL cp_add_iter_level(logger%iter_info,level_name="EP_LIN_SOLVER",error=error) 00890 CALL cp_iterate(logger%iter_info,error=error) 00891 DO i=1,20 00892 ! orthogonalize x wrt. to psi0 00893 CALL p_preortho(p_env=ep_env%main_p_env, & 00894 qs_env=ep_env%main_qs_env, v=ep_env%psi1,& 00895 n_cols=ep_env%main_p_env%n_mo, error=error) 00896 ! localize x 00897 CALL ep_env_localize_matrix(ep_env, ep_env%psi1, error=error) 00898 ! r=b-A x 00899 CALL p_op_ep(p_env=ep_env%main_p_env, qs_env=ep_env%main_qs_env,& 00900 v=ep_env%psi1, res=r, error=error) 00901 DO ispin=1,nspins 00902 CALL cp_fm_scale_and_add(alpha=-1.0_dp,matrix_a=r(ispin)%matrix,& 00903 beta=1.0_dp,matrix_b=ep_env%m_pi_Hrho_psi0d(ispin)%matrix,& 00904 error=error) 00905 END DO 00906 ! localize r 00907 CALL ep_env_localize_matrix(ep_env, r, error=error) 00908 00909 ! check convergence (with ||r||) 00910 r_norm2=0.0_dp 00911 DO ispin=1,nspins 00912 CALL cp_fm_trace(matrix_a=r(ispin)%matrix,& 00913 matrix_b=r(ispin)%matrix,trace=r_norm2_s,error=error) 00914 r_norm2=r_norm2+r_norm2_s 00915 END DO 00916 IF (r_norm2<eps_r**2) EXIT 00917 00918 ! z=M^-1 r 00919 DO ispin=1,nspins 00920 CALL apply_preconditioner(ep_env%precond(ispin)%preconditioner,& 00921 matrix_in=r(ispin)%matrix,& 00922 matrix_out=z(ispin)%matrix,error=error) 00923 END DO 00924 ! localize z 00925 CALL ep_env_localize_matrix(ep_env, z, error=error) 00926 00927 ! p=z 00928 DO ispin=1,nspins 00929 CALL cp_fm_scale_and_add(alpha=0.0_dp,matrix_a=p(ispin)%matrix,beta=1.0_dp,& 00930 matrix_b=z(ispin)%matrix,error=error) 00931 END DO 00932 r_z=0.0_dp 00933 DO ispin=1,nspins 00934 CALL cp_fm_trace(matrix_a=r(ispin)%matrix,& 00935 matrix_b=z(ispin)%matrix,trace=r_z_s,error=error) 00936 r_z=r_z+r_z_s 00937 END DO 00938 00939 DO j=1,10 00940 CALL cp_iterate(logger%iter_info,error=error) 00941 CALL cp_log(logger,cp_note_level,routineP,"i="//cp_to_string(i)//& 00942 "j="//cp_to_string(j)//"r_norm2="//cp_to_string(r_norm2)) 00943 00944 ! put the component of p orthogonal to psi0 in p_ortho 00945 DO ispin=1,nspins 00946 CALL cp_fm_scale_and_add(alpha=0.0_dp, matrix_a=p_ortho(ispin)%matrix,& 00947 beta=1.0_dp,& 00948 matrix_b=p(ispin)%matrix,error=error) 00949 END DO 00950 CALL p_preortho(p_env=ep_env%main_p_env,& 00951 qs_env=ep_env%main_qs_env, v=p_ortho,& 00952 n_cols=ep_env%main_p_env%n_mo, error=error) 00953 ! localize p_ortho 00954 CALL ep_env_localize_matrix(ep_env, p_ortho, error=error) 00955 00956 ! calc norm of p_ortho 00957 p_ortho_norm=0.0_dp 00958 DO ispin=1,nspins 00959 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix,& 00960 p_ortho(ispin)%matrix,Ap(ispin)%matrix,& 00961 ncol=ep_env%main_p_env%n_mo(ispin),& 00962 error=error) 00963 CALL cp_fm_trace(matrix_a=p_ortho(ispin)%matrix,& 00964 matrix_b=Ap(ispin)%matrix,& 00965 trace=p_ortho_norm_s,error=error) 00966 p_ortho_norm=p_ortho_norm+p_ortho_norm_s 00967 END DO 00968 p_ortho_norm=SQRT(p_ortho_norm) 00969 00970 ! rescale p_ortho 00971 CPPreconditionNoFail(p_ortho_norm>eps_r*0.5,cp_warning_level,routineP,error) 00972 DO ispin=1,nspins 00973 CALL cp_fm_scale_and_add(alpha=1.0_dp/p_ortho_norm,& 00974 matrix_a=p_ortho(ispin)%matrix,error=error) 00975 END DO 00976 00977 ! Ap=A p_ortho=p_ortho_norm A p 00978 CALL p_op_ep(p_env=ep_env%main_p_env, qs_env=ep_env%main_qs_env,& 00979 v=p_ortho, res=Ap,error=error) 00980 ! localize Ap 00981 CALL ep_env_localize_matrix(ep_env, Ap, error=error) 00982 00983 Ap_p=0.0_dp 00984 DO ispin=1,nspins 00985 CALL cp_fm_trace(matrix_a=Ap(ispin)%matrix,& 00986 matrix_b=p(ispin)%matrix,trace=Ap_p_s,error=error) 00987 Ap_p=Ap_p+Ap_p_s 00988 END DO 00989 00990 ! this alpha is the alpha in alg. desc. times p_ortho_norm 00991 alpha=r_z/Ap_p 00992 00993 ! x=x+alpha p_ortho 00994 DO ispin=1,nspins 00995 CALL cp_fm_scale_and_add(alpha=1.0_dp,matrix_a=ep_env%psi1(ispin)%matrix,& 00996 beta=alpha,matrix_b=p_ortho(ispin)%matrix,& 00997 error=error) 00998 END DO 00999 01000 ! r=r- alpha Ap 01001 DO ispin=1,nspins 01002 CALL cp_fm_scale_and_add(alpha=1.0_dp,matrix_a=r(ispin)%matrix,& 01003 beta=-alpha,matrix_b=Ap(ispin)%matrix,& 01004 error=error) 01005 END DO 01006 01007 r_norm2=0.0_dp 01008 DO ispin=1,nspins 01009 CALL cp_fm_trace(matrix_a=r(ispin)%matrix,& 01010 matrix_b=r(ispin)%matrix,trace=r_norm2_s,error=error) 01011 r_norm2=r_norm2+r_norm2_s 01012 END DO 01013 IF (r_norm2<eps_r**2) EXIT 01014 01015 ! z=M^-1 r 01016 DO ispin=1,nspins 01017 CALL apply_preconditioner(ep_env%precond(ispin)%preconditioner,& 01018 matrix_in=r(ispin)%matrix,& 01019 matrix_out=z(ispin)%matrix,error=error) 01020 END DO 01021 ! localize z 01022 CALL ep_env_localize_matrix(ep_env, z, error=error) 01023 01024 r_z_new=0.0_dp 01025 DO ispin=1,nspins 01026 CALL cp_fm_trace(matrix_a=r(ispin)%matrix,& 01027 matrix_b=z(ispin)%matrix,trace=r_z_s,error=error) 01028 r_z_new=r_z_new+r_z_s 01029 END DO 01030 01031 beta= r_z_new/r_z 01032 r_z=r_z_new 01033 01034 ! p=beta p+z 01035 DO ispin=1,nspins 01036 CALL cp_fm_scale_and_add(alpha=beta,matrix_a=p(ispin)%matrix,& 01037 beta=1.0_dp,matrix_b=z(ispin)%matrix,& 01038 error=error) 01039 END DO 01040 END DO 01041 01042 output_unit=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%RUN_INFO/LIN_SOLV",& 01043 extension=".epLog",error=error) 01044 IF (output_unit>0) THEN 01045 WRITE(output_unit,"(a,a,i6,a,f20.15)")routineP,"i=",i,& 01046 "j=final r_norm2=",r_norm2 01047 CALL cp_print_key_finished_output(output_unit,logger,ep_env%input,& 01048 "PRINT%RUN_INFO/LIN_SOLV",& 01049 error=error) 01050 END IF 01051 END DO 01052 01053 CALL cp_rm_iter_level(logger%iter_info,level_name="EP_LIN_SOLVER",error=error) 01054 output_unit=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%RUN_INFO/LIN_SOLV",& 01055 extension=".epLog",error=error) 01056 IF (output_unit>0) THEN 01057 WRITE(output_unit,"(a,f20.15)") "solve finished, r_norm2=",r_norm2 01058 CALL cp_print_key_finished_output(output_unit,logger,ep_env%input,& 01059 "PRINT%RUN_INFO/LIN_SOLV",& 01060 error=error) 01061 END IF 01062 CPPostcondition(r_norm2<eps_r**2,cp_warning_level,routineP,error,failure) 01063 01064 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,p,error=error) 01065 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,Ap,error=error) 01066 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,r,error=error) 01067 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,z,error=error) 01068 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,p_ortho,error=error) 01069 END IF 01070 01071 END SUBROUTINE stupid_solve 01072 01073 ! ***************************************************************************** 01086 SUBROUTINE p_op_ep(p_env, qs_env, v, res, error) 01087 TYPE(qs_p_env_type), POINTER :: p_env 01088 TYPE(qs_environment_type), POINTER :: qs_env 01089 TYPE(cp_fm_p_type), DIMENSION(:), 01090 INTENT(in) :: v 01091 TYPE(cp_fm_p_type), DIMENSION(:), 01092 INTENT(inout) :: res 01093 TYPE(cp_error_type), INTENT(inout) :: error 01094 01095 CHARACTER(len=*), PARAMETER :: routineN = 'p_op_ep', 01096 routineP = moduleN//':'//routineN 01097 01098 INTEGER :: handle, ispin, lfomo, nmo 01099 LOGICAL :: failure 01100 REAL(KIND=dp) :: maxocc, scale_factor 01101 TYPE(mo_set_p_type), DIMENSION(:), 01102 POINTER :: mos 01103 01104 CALL timeset(routineN,handle) 01105 01106 failure=.FALSE. 01107 NULLIFY(mos) 01108 01109 CPPrecondition(ASSOCIATED(p_env),cp_failure_level,routineP,error,failure) 01110 CPPrecondition(p_env%ref_count>0,cp_failure_level,routineP,error,failure) 01111 IF (.NOT. failure) THEN 01112 CALL get_qs_env(qs_env,mos=mos,error=error) 01113 CALL p_op_l1(p_env=p_env, qs_env=qs_env,& 01114 v=v, res=res,error=error) 01115 CALL get_mo_set(mos(1)%mo_set,lfomo=lfomo, nmo=nmo, maxocc=maxocc) 01116 IF (lfomo>nmo) THEN 01117 scale_factor=maxocc 01118 DO ispin=2,SIZE(mos) 01119 CALL get_mo_set(mos(ispin)%mo_set,lfomo=lfomo, & 01120 nmo=nmo, maxocc=maxocc) 01121 CPAssert(scale_factor==maxocc,cp_failure_level,routineP,error,failure) 01122 END DO 01123 CALL p_op_l2_fawzi(p_env=p_env, qs_env=qs_env,& 01124 v=v, res=res,& 01125 alpha=scale_factor,beta=scale_factor,error=error) 01126 ELSE 01127 CALL cp_unimplemented_error(fromWhere=routineP,& 01128 message="symmetrized onesided smearing to do",& 01129 error=error) 01130 END IF 01131 CALL p_postortho(p_env=p_env, qs_env=qs_env, v=res, n_cols=p_env%n_mo,& 01132 error=error) 01133 END IF 01134 01135 CALL timestop(handle) 01136 01137 END SUBROUTINE p_op_ep 01138 01139 ! ***************************************************************************** 01151 SUBROUTINE p_env_write_ep_matrix(p_env, qs_env, name, error) 01152 TYPE(qs_p_env_type), POINTER :: p_env 01153 TYPE(qs_environment_type), POINTER :: qs_env 01154 CHARACTER(len=*), INTENT(in) :: name 01155 TYPE(cp_error_type), INTENT(inout) :: error 01156 01157 CHARACTER(len=*), PARAMETER :: routineN = 'p_env_write_ep_matrix', 01158 routineP = moduleN//':'//routineN 01159 01160 INTEGER :: iao, imo, ispin, ispin2 01161 LOGICAL :: failure 01162 REAL(KIND=dp) :: v_norm, v_norm_s 01163 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01164 POINTER :: matrix_s 01165 TYPE(cp_fm_p_type), DIMENSION(:), 01166 POINTER :: res, v 01167 TYPE(cp_fm_pool_p_type), DIMENSION(:), 01168 POINTER :: ao_mo_fm_pools 01169 TYPE(cp_logger_type), POINTER :: logger 01170 TYPE(dft_control_type), POINTER :: dft_control 01171 01172 failure=.FALSE. 01173 logger => cp_error_get_logger(error) 01174 NULLIFY(ao_mo_fm_pools, dft_control,matrix_s) 01175 01176 CPPrecondition(ASSOCIATED(p_env),cp_failure_level,routineP,error,failure) 01177 CPPrecondition(p_env%ref_count>0,cp_failure_level,routineP,error,failure) 01178 IF (.NOT. failure) THEN 01179 CALL get_qs_env(qs_env,& 01180 dft_control=dft_control,matrix_s=matrix_s,error=error) 01181 01182 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools,& 01183 error=error) 01184 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,v,error=error) 01185 CALL fm_pools_create_fm_vect(ao_mo_fm_pools,res,name=name,& 01186 error=error) 01187 01188 DO ispin=1,dft_control%nspins 01189 DO imo=1,p_env%n_mo(ispin) 01190 DO iao=1,p_env%n_ao(ispin) 01191 CALL cp_fm_vect_set_all(v,0.0_dp,error=error) 01192 CALL cp_fm_set_element(v(ispin)%matrix,iao,imo,1.0_dp,& 01193 error=error) 01194 01195 CALL p_preortho(p_env=p_env, qs_env=qs_env, v=v,& 01196 n_cols=p_env%n_mo, error=error) 01197 ! scale v 01198 v_norm=0.0_dp 01199 DO ispin2=1,dft_control%nspins 01200 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix,& 01201 v(ispin2)%matrix,res(ispin2)%matrix,& 01202 ncol=p_env%n_mo(ispin2), & 01203 error=error) 01204 CALL cp_fm_trace(matrix_a=v(ispin2)%matrix, & 01205 matrix_b=res(ispin2)%matrix, trace=v_norm_s, error=error) 01206 v_norm=v_norm+v_norm_s 01207 END DO 01208 v_norm=SQRT(v_norm) 01209 CPPreconditionNoFail(v_norm>1.0e-8,cp_warning_level,routineP,error) 01210 DO ispin2=1,dft_control%nspins 01211 CALL cp_fm_scale_and_add(1.0_dp/v_norm,v(ispin2)%matrix,error=error) 01212 END DO 01213 01214 CALL p_op_ep(p_env=p_env, qs_env=qs_env, v=v,& 01215 res=res, error=error) 01216 01217 DO ispin2=1,dft_control%nspins 01218 CALL cp_fm_scale_and_add(v_norm,res(ispin2)%matrix,error=error) 01219 END DO 01220 01221 CALL cp_fm_vect_write(matrixes=res, & 01222 unit_nr=cp_logger_get_default_unit_nr(logger),& 01223 long_description=.TRUE., error=error) 01224 CALL m_flush(cp_logger_get_default_unit_nr(logger)) 01225 END DO 01226 END DO 01227 END DO 01228 01229 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools, v, error=error) 01230 CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools, res, error=error) 01231 END IF 01232 01233 01234 END SUBROUTINE p_env_write_ep_matrix 01235 01236 ! ***************************************************************************** 01244 SUBROUTINE ep_env_localize_matrix(ep_env, wf_coeff, error) 01245 TYPE(ep_env_type), POINTER :: ep_env 01246 TYPE(cp_fm_p_type), DIMENSION(:), 01247 POINTER :: wf_coeff 01248 TYPE(cp_error_type), INTENT(inout) :: error 01249 01250 CHARACTER(len=*), PARAMETER :: routineN = 'ep_env_localize_matrix', 01251 routineP = moduleN//':'//routineN 01252 01253 INTEGER :: handle, ispin, isub, nmo, 01254 nspins, start_col 01255 LOGICAL :: failure 01256 TYPE(cp_logger_type), POINTER :: logger 01257 TYPE(qs_environment_type), POINTER :: main_qs_env 01258 TYPE(qs_p_projection_p_type), 01259 DIMENSION(:), POINTER :: sub_proj 01260 01261 CALL timeset(routineN,handle) 01262 01263 failure=.FALSE. 01264 NULLIFY(main_qs_env, sub_proj) 01265 logger => cp_error_get_logger(error) 01266 01267 CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 01268 CALL ep_env_get(ep_env,& 01269 sub_proj=sub_proj,main_qs_env=main_qs_env,error=error) 01270 CALL cp_error_check(error,failure) 01271 CPPrecondition(ASSOCIATED(sub_proj),cp_failure_level,routineP,error,failure) 01272 CPPrecondition(ASSOCIATED(main_qs_env),cp_failure_level,routineP,error,failure) 01273 IF (.NOT. failure) THEN 01274 nspins=main_qs_env%dft_control%nspins 01275 DO ispin=1,nspins 01276 start_col=1 01277 CALL cp_fm_get_info(wf_coeff(ispin)%matrix,ncol_global=nmo,error=error) 01278 01279 DO isub=1,SIZE(sub_proj) 01280 CPInvariant(start_col+ep_env%sub_nmo(ispin)<nmo+2,cp_failure_level,routineP,error,failure) 01281 start_col=start_col+ep_env%sub_nmo(ispin) 01282 END DO 01283 CPPostcondition(start_col==nmo+1,cp_failure_level,routineP,error,failure) 01284 END DO 01285 END IF 01286 01287 CALL timestop(handle) 01288 END SUBROUTINE ep_env_localize_matrix 01289 01290 ! ***************************************************************************** 01297 SUBROUTINE ep_env_calc_e_f_low(ep_env_id, calc_f, ierr) 01298 INTEGER, INTENT(in) :: ep_env_id, calc_f 01299 INTEGER, INTENT(out) :: ierr 01300 01301 CHARACTER(len=*), PARAMETER :: routineN = 'ep_env_calc_e_f_low', 01302 routineP = moduleN//':'//routineN 01303 01304 LOGICAL :: failure 01305 TYPE(cp_error_type) :: error 01306 TYPE(ep_env_type), POINTER :: ep_env 01307 01308 failure=.FALSE. 01309 01310 CALL cp_error_init(error) 01311 ep_env => ep_envs_get_ep_env(ep_env_id) 01312 CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 01313 IF (.NOT.failure) THEN 01314 CALL ep_env_calc_e0(ep_env,error=error) 01315 END IF 01316 CALL cp_error_check(error,failure) 01317 IF (failure) THEN 01318 ierr=1 01319 ELSE 01320 ierr=0 01321 END IF 01322 CALL cp_error_dealloc_ref(error) 01323 END SUBROUTINE ep_env_calc_e_f_low 01324 01325 ! ***************************************************************************** 01335 SUBROUTINE qs_check_i_alloc(qs_env,error) 01336 TYPE(qs_environment_type), POINTER :: qs_env 01337 TYPE(cp_error_type), INTENT(inout) :: error 01338 01339 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_check_i_alloc', 01340 routineP = moduleN//':'//routineN 01341 01342 INTEGER :: handle, homo, ispin 01343 LOGICAL :: failure, 01344 gth_potential_present, 01345 uniform_occupation 01346 TYPE(atomic_kind_type), DIMENSION(:), 01347 POINTER :: atomic_kind_set 01348 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01349 POINTER :: matrix_ks, matrix_s 01350 TYPE(cp_fm_pool_p_type), DIMENSION(:), 01351 POINTER :: ao_mo_fm_pools 01352 TYPE(cp_fm_type), POINTER :: mo_coeff 01353 TYPE(cp_logger_type), POINTER :: logger 01354 TYPE(dft_control_type), POINTER :: dft_control 01355 TYPE(mo_set_p_type), DIMENSION(:), 01356 POINTER :: mos 01357 TYPE(qs_ks_env_type), POINTER :: ks_env 01358 01359 ! my_transition_potential, & 01360 01361 CALL timeset(routineN,handle) 01362 01363 NULLIFY(matrix_ks, ao_mo_fm_pools, matrix_s, dft_control, mos, & 01364 ks_env) 01365 NULLIFY(atomic_kind_set, mo_coeff) 01366 logger => cp_error_get_logger(error) 01367 failure=.FALSE. 01368 01369 ! my_transition_potential = .FALSE. 01370 uniform_occupation = .TRUE. 01371 01372 CALL get_qs_env(qs_env=qs_env,& 01373 dft_control=dft_control,& 01374 mos=mos,& 01375 matrix_ks=matrix_ks,& 01376 ks_env=ks_env,& 01377 atomic_kind_set=atomic_kind_set,& 01378 matrix_s=matrix_s,error=error) 01379 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01380 gth_potential_present=gth_potential_present) 01381 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools,& 01382 error=error) 01383 01384 01385 ! *** finish initialization of the MOs *** 01386 CPPrecondition(ASSOCIATED(mos),cp_failure_level,routineP,error,failure) 01387 IF (.NOT.failure) THEN 01388 DO ispin=1,SIZE(mos) 01389 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff,homo=homo) 01390 IF (.NOT.ASSOCIATED(mo_coeff)) THEN 01391 CALL init_mo_set(mos(ispin)%mo_set,& 01392 ao_mo_fm_pools(ispin)%pool,& 01393 name="qs_env"//TRIM(ADJUSTL(cp_to_string(qs_env%id_nr)))//& 01394 "%mo"//TRIM(ADJUSTL(cp_to_string(ispin))),& 01395 error=error) 01396 END IF 01397 ! IF(my_transition_potential .AND. ispin==1) THEN 01398 ! CALL get_mo_set(mos(ispin)%mo_set,& 01399 ! occupation_numbers=occupation_numbers) 01400 ! occupation_numbers(dft_control%xas_estate) = & 01401 ! dft_control%xas_control%occ_estate 01402 ! occupation_numbers(homo) = dft_control%xas_control%occ_homo 01403 ! uniform_occupation = .FALSE. 01404 ! 01405 ! CALL set_mo_set(mos(ispin)%mo_set,& 01406 ! uniform_occupation=uniform_occupation,error=error) 01407 ! END IF 01408 END DO 01409 END IF 01410 01411 ! *** get the mo_derivs OK if needed *** 01412 IF (qs_env%requires_mo_derivs) THEN 01413 ! this is never used ... 01414 !CALL get_qs_env(qs_env,mo_derivs=mo_derivs,error=error) 01415 !IF (.NOT.ASSOCIATED(mo_derivs)) THEN 01416 ! IF (dft_control%restricted) THEN ! right now, there might be more mos than needed derivs 01417 ! ALLOCATE(mo_derivs(1)) 01418 ! CALL get_mo_set(mos(1)%mo_set,mo_coeff=mo_coeff) 01419 ! CALL cp_fm_create(mo_derivs(1)%matrix,mo_coeff%matrix_struct,error=error) 01420 ! ELSE 01421 ! ALLOCATE(mo_derivs(dft_control%nspins)) 01422 ! DO ispin=1,dft_control%nspins 01423 ! CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff) 01424 ! CALL cp_fm_create(mo_derivs(ispin)%matrix,mo_coeff%matrix_struct,error=error) 01425 ! ENDDO 01426 ! ENDIF 01427 ! CALL set_qs_env(qs_env,mo_derivs=mo_derivs,error=error) 01428 !ENDIF 01429 ELSE 01430 ! nothing should be done 01431 ENDIF 01432 01433 ! *** Allocate matrix_ks and put it in the QS environment *** 01434 01435 IF (.not.ASSOCIATED(matrix_ks)) THEN 01436 CALL cp_dbcsr_allocate_matrix_set( matrix_ks, dft_control%nspins, error ) 01437 DO ispin=1,SIZE(matrix_ks) 01438 ALLOCATE(matrix_ks(ispin)%matrix) 01439 CALL cp_dbcsr_init(matrix_ks(ispin)%matrix,error=error) 01440 CALL cp_dbcsr_copy(matrix_ks(ispin)%matrix, matrix_s(ispin)%matrix,& 01441 name="EP"//TRIM(cp_iter_string(logger%iter_info,error=error))//& 01442 "KOHN-SHAM_MATRIX", error=error) 01443 ENDDO 01444 01445 CALL set_qs_env(qs_env=qs_env,& 01446 matrix_ks=matrix_ks,error=error) 01447 01448 END IF 01449 01450 ! *** allocate the ks env ** 01451 IF (.NOT.ASSOCIATED(ks_env)) THEN 01452 CALL qs_ks_create(ks_env,qs_env=qs_env,error=error) 01453 CALL set_qs_env(qs_env, ks_env=ks_env,error=error) 01454 CALL qs_ks_release(ks_env,error=error) 01455 END IF 01456 01457 CALL timestop(handle) 01458 01459 END SUBROUTINE qs_check_i_alloc 01460 01461 ! ***************************************************************************** 01466 SUBROUTINE ep_env_calc_e0(ep_env,error) 01467 TYPE(ep_env_type), POINTER :: ep_env 01468 TYPE(cp_error_type), INTENT(inout) :: error 01469 01470 CHARACTER(len=*), PARAMETER :: routineN = 'ep_env_calc_e0', 01471 routineP = moduleN//':'//routineN 01472 01473 INTEGER :: flag, handle, i, iat, iatom, idir, ii, ikind, imol, ispin, iw, 01474 lfomo, natom, nkind, nmo, nspin, output_unit, stat, unit_nr 01475 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 01476 natom_of_kind 01477 INTEGER, DIMENSION(:), POINTER :: at 01478 LOGICAL :: do_rotations, e0_only, failure 01479 REAL(dp) :: e_corr_s, eps_lin_solv, maxocc 01480 TYPE(atomic_kind_type), DIMENSION(:), 01481 POINTER :: atomic_kind_set 01482 TYPE(cp_2d_r_p_type), ALLOCATABLE, 01483 DIMENSION(:) :: rotatedC0 01484 TYPE(cp_blacs_env_type), POINTER :: blacs_env 01485 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01486 POINTER :: matrix_ks, matrix_s, matrix_w 01487 TYPE(cp_error_type) :: suberror 01488 TYPE(cp_fm_p_type), ALLOCATABLE, 01489 DIMENSION(:) :: psi0 01490 TYPE(cp_fm_type), POINTER :: orbitals 01491 TYPE(cp_logger_type), POINTER :: logger 01492 TYPE(cp_para_env_type), POINTER :: para_env 01493 TYPE(cp_subsys_type), POINTER :: subsys 01494 TYPE(dft_control_type), POINTER :: dft_control 01495 TYPE(global_environment_type), POINTER :: globenv 01496 TYPE(mo_set_p_type), DIMENSION(:), 01497 POINTER :: mos 01498 TYPE(particle_list_type), POINTER :: particles 01499 TYPE(particle_type), DIMENSION(:), 01500 POINTER :: particle_set 01501 TYPE(qs_energy_type), POINTER :: energy 01502 TYPE(qs_environment_type), POINTER :: qs_env 01503 TYPE(qs_force_type), DIMENSION(:), 01504 POINTER :: force 01505 TYPE(qs_ks_env_type), POINTER :: ks_env 01506 TYPE(qs_rho_type), POINTER :: rho 01507 TYPE(replica_env_type), POINTER :: rep_env 01508 TYPE(scf_control_type), POINTER :: scf_control 01509 TYPE(section_vals_type), POINTER :: print_section, qs_env_input 01510 01511 CALL timeset(routineN,handle) 01512 failure = .FALSE. 01513 NULLIFY (logger) 01514 logger => cp_error_get_logger(error) 01515 01516 NULLIFY (atomic_kind_set, dft_control, force, ks_env, mos, matrix_ks, & 01517 matrix_s, matrix_w, particle_set, rho, scf_control,para_env,& 01518 blacs_env,orbitals) 01519 01520 globenv => ep_env%globenv 01521 qs_env => ep_env%main_qs_env 01522 para_env=>qs_env%para_env 01523 01524 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01525 extension=".epLog",message=routineP//& 01526 " started",error=error) 01527 CALL cp_error_init(suberror,template_error=error) 01528 01529 CALL get_qs_env(qs_env=qs_env,& 01530 atomic_kind_set=atomic_kind_set,& 01531 dft_control=dft_control,& 01532 force=force,& 01533 particle_set=particle_set,& 01534 scf_control=scf_control,& 01535 ks_env=ks_env,& 01536 input=qs_env_input,& 01537 energy=energy,error=error) 01538 CALL get_qs_env(qs_env=qs_env,subsys=subsys,error=error) 01539 CALL cp_subsys_get(subsys,particles=particles,error=error) 01540 01541 01542 ! ensure storage of electronic kinetic energy 01543 CALL section_vals_val_get(qs_env_input,"DFT%PRINT%AO_MATRICES%__CONTROL_VAL",& 01544 i_val=flag,error=error) 01545 flag=IBSET(flag,cp_p_store) 01546 CALL section_vals_val_set(qs_env_input,"DFT%PRINT%AO_MATRICES%__CONTROL_VAL",& 01547 i_val=flag,error=error) 01548 CALL section_vals_val_set(qs_env_input,"DFT%PRINT%AO_MATRICES%KINETIC_ENERGY",& 01549 l_val=.TRUE.,error=error) 01550 01551 natom = SIZE(particle_set) 01552 01553 ! zero out the forces 01554 DO iatom=1,natom 01555 particle_set(iatom)%f=0.0_dp 01556 END DO 01557 ! zero energy 01558 CALL ep_energy_zero(ep_env%energy,error=error) 01559 01560 ALLOCATE (atom_of_kind(natom),STAT=stat) 01561 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01562 01563 ALLOCATE (kind_of(natom),STAT=stat) 01564 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01565 01566 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01567 atom_of_kind=atom_of_kind,& 01568 kind_of=kind_of) 01569 01570 IF (.NOT.ASSOCIATED(force)) THEN 01571 ! *** Allocate the force data structure *** 01572 nkind = SIZE(atomic_kind_set) 01573 ALLOCATE (natom_of_kind(nkind),STAT=stat) 01574 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01575 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01576 natom_of_kind=natom_of_kind) 01577 CALL allocate_qs_force(force,natom_of_kind) 01578 DEALLOCATE (natom_of_kind,STAT=stat) 01579 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01580 CALL set_qs_env(qs_env,force=force,error=error) 01581 END IF 01582 CALL zero_qs_force(force) 01583 CALL ep_force_zero(ep_env%force,error=error) 01584 01585 ! !!!!! qs energy 01586 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01587 extension=".epLog",message=routineP//& 01588 " qs_energy setup in main_qs_env",error=error) 01589 CALL build_qs_neighbor_lists(qs_env,para_env,force_env_section=qs_env%input,error=error) 01590 01591 ! *** Calculate the overlap and the core Hamiltonian integral matrix *** 01592 CALL build_core_hamiltonian_matrix(qs_env=qs_env,calculate_forces=.FALSE.,error=error) 01593 CALL qs_env_update_s_mstruct(qs_env,error=error) 01594 CALL qs_check_i_alloc(qs_env,error=error) 01595 CALL calculate_ecore_self(qs_env,error=error) 01596 CALL calculate_ecore_overlap(qs_env, para_env, & 01597 calculate_forces=.FALSE.,error=error) 01598 01599 ! *** calculate psi0 01600 ! ** update positions of replicas 01601 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01602 extension=".epLog",message=routineP//& 01603 " calculating psi0",error=error) 01604 rep_env => ep_env%mol_envs 01605 CALL get_qs_env(ep_env%main_qs_env,particle_set=particle_set,& 01606 para_env=para_env,blacs_env=blacs_env,mos=mos,error=error) 01607 DO imol=1,ep_env%nmol 01608 at => ep_env%sub_proj(imol)%projection%atoms 01609 ii=0 01610 DO iat=1,SIZE(at) 01611 DO idir=1,3 01612 ii=ii+1 01613 rep_env%r(ii,imol)=particle_set(at(iat))%r(idir) 01614 END DO 01615 END DO 01616 END DO 01617 01618 CALL section_vals_val_get(ep_env%root_section,"FORCE_EVAL%EP%ROTATE",& 01619 l_val=do_rotations,error=error) 01620 IF (.NOT.do_rotations) THEN 01621 ! ** calc e and psi00 01622 IF (ASSOCIATED(rep_env)) THEN 01623 CALL rep_env_calc_e_f(rep_env,calc_f=.FALSE.,error=error) 01624 END IF 01625 ep_env%energy%e_no_int=SUM(rep_env%f(rep_env%ndim+1,:)) 01626 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01627 extension=".epLog",message=routineP//& 01628 " calculated psi0 locally, e_no_int="//& 01629 TRIM(ADJUSTL(cp_to_string(ep_env%energy%e_no_int))),error=error) 01630 ! ** transfer psi0 01631 CALL ep_env_transfer_psi0(ep_env,error=error) 01632 ELSE 01633 ! matrix for rotated C 01634 ALLOCATE(rotatedC0(ep_env%nspins),stat=stat) 01635 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01636 DO ispin=1,ep_env%nspins 01637 ALLOCATE(rotatedC0(ispin)%array(ep_env%sub_nmo(ispin),ep_env%sub_nao(ispin)),stat=stat) 01638 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01639 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=orbitals,nmo=nmo) 01640 CALL cp_fm_set_all(orbitals,0.0_dp,error=error) 01641 END DO 01642 DO imol=1,ep_env%nmol 01643 IF (rep_env%rep_is_local(imol)) THEN 01644 ! this should return the C0 for a molecule 01645 CALL calc_c0_tilde(ep_env,& 01646 mol_pos=rep_env%r(:,imol),new_c=rotatedC0,& 01647 imol=imol,error=error) 01648 END IF 01649 DO ispin=1,ep_env%nspins 01650 ! transfer to matrix in the full systems (psi0) 01651 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=orbitals,nmo=nmo) 01652 CALL mp_bcast(rotatedC0(ispin)%array, rep_env%replica_owner(imol),& 01653 rep_env%para_env_inter_rep%group) 01654 CALL cp_fm_set_submatrix(orbitals,rotatedC0(ispin)%array,& 01655 start_col=(imol-1)*ep_env%sub_nmo(ispin)+1,& 01656 start_row=(imol-1)*ep_env%sub_nao(ispin)+1,& 01657 transpose=.TRUE.,error=error) 01658 END DO 01659 END DO 01660 DO ispin=1,ep_env%nspins 01661 DEALLOCATE(rotatedC0(ispin)%array,stat=stat) 01662 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 01663 END DO 01664 DEALLOCATE(rotatedC0,stat=stat) 01665 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 01666 END IF 01667 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01668 extension=".epLog",message=routineP//& 01669 " transferred psi0 to main system",error=error) 01670 01671 ! update main_p_env & psi0d 01672 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01673 extension=".epLog",message=routineP//& 01674 " update main_p_env",error=error) 01675 CALL p_env_did_change(ep_env%sub_p_env, s_struct_changed=.TRUE., & 01676 grid_changed=.FALSE.,error=error) 01677 CALL p_env_psi0_changed(ep_env%main_p_env, qs_env=ep_env%main_qs_env,& 01678 Hrho_psi0d=ep_env%m_pi_Hrho_psi0d, error=error) 01679 CALL p_postortho(p_env=ep_env%main_p_env,qs_env=ep_env%main_qs_env,& 01680 v=ep_env%m_pi_Hrho_psi0d, & 01681 n_cols=ep_env%main_p_env%n_mo,& 01682 error=error) 01683 01684 ! ** update qs_env%rho 01685 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01686 extension=".epLog",message=routineP//& 01687 " update qs_env%rho and ks matrix",error=error) 01688 CALL get_qs_env(qs_env, rho=rho, mos=mos, error=error) 01689 DO ispin=1,ep_env%nspins 01690 CALL cp_dbcsr_set(rho%rho_ao(ispin)%matrix,0._dp,error=error) 01691 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho%rho_ao(ispin)%matrix,& 01692 matrix_v=mos(ispin)%mo_set%mo_coeff,& 01693 matrix_g=ep_env%main_p_env%psi0d(ispin)%matrix,& 01694 ncol=ep_env%full_nmo(ispin),& 01695 error=error) 01696 CALL get_mo_set(mos(ispin)%mo_set, maxocc=maxocc) 01697 CALL cp_dbcsr_scale(rho%rho_ao(ispin)%matrix,alpha_scalar=maxocc,error=error) 01698 END DO 01699 CALL qs_rho_update_rho(rho, qs_env=qs_env, error=error) 01700 CALL qs_ks_did_change(ks_env,rho_changed=.TRUE.,error=error) 01701 CALL qs_ks_update_qs_env(ks_env,qs_env,calculate_forces=.FALSE.,& 01702 just_energy=.FALSE.,error=error) 01703 ep_env%energy%e0=ep_env%main_qs_env%energy%total 01704 unit_nr=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%RUN_INFO",& 01705 extension=".epLog",error=error) 01706 CALL print_qs_energies(qs_env,unit_nr,error) 01707 IF (unit_nr>0) CALL cp_print_key_finished_output(unit_nr,& 01708 logger,ep_env%input,"PRINT%RUN_INFO", error=error) 01709 01710 DO ispin=1,SIZE(mos) 01711 CALL get_mo_set(mos(ispin)%mo_set,lfomo=lfomo, & 01712 nmo=nmo, maxocc=maxocc) 01713 IF (lfomo>nmo) THEN 01714 CALL cp_fm_scale_and_add(alpha=-maxocc,& 01715 matrix_a=ep_env%m_pi_Hrho_psi0d(ispin)%matrix,error=error) 01716 ELSE 01717 CALL cp_unimplemented_error(fromWhere=routineP,& 01718 message="symmetrized onesided smearing to do",& 01719 error=error) 01720 END IF 01721 END DO 01722 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01723 extension=".epLog",message=routineP//& 01724 " finished qs_enegy setup",error=error) 01725 ! !!!! end qs_energy 01726 01727 ! init preconditioner 01728 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01729 extension=".epLog",message=routineP//& 01730 " init proeconditioner main_p_env",error=error) 01731 ALLOCATE(psi0(dft_control%nspins),stat=stat) 01732 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01733 IF (.NOT.failure) THEN 01734 DO ispin=1,dft_control%nspins 01735 NULLIFY(psi0(ispin)%matrix) 01736 CALL get_mo_set(ep_env%main_qs_env%mos(ispin)%mo_set,& 01737 mo_coeff=psi0(ispin)%matrix) 01738 END DO 01739 ALLOCATE(ep_env%precond(dft_control%nspins),stat=stat) 01740 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01741 ispin=1 01742 NULLIFY(ep_env%precond(ispin)%preconditioner) 01743 ALLOCATE(ep_env%precond(ispin)%preconditioner,stat=stat) 01744 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01745 CALL init_preconditioner(ep_env%precond(ispin)%preconditioner,& 01746 para_env=para_env,& 01747 blacs_env=blacs_env,error=error) 01748 CALL make_preconditioner(ep_env%precond(ispin)%preconditioner,& 01749 ot_precond_full_kinetic,& 01750 ot_precond_solver_default,& 01751 ep_env%main_qs_env%matrix_ks(ispin)%matrix,& 01752 ep_env%main_qs_env%matrix_s(1)%matrix,& 01753 ep_env%main_qs_env%kinetic(1)%matrix,& 01754 ep_env%main_qs_env%mos(ispin)%mo_set,0.2_dp,error=error) 01755 ! same precond for both spins 01756 DO ispin=2,dft_control%nspins 01757 ep_env%precond(ispin)%preconditioner => & 01758 ep_env%precond(1)%preconditioner 01759 END DO 01760 END IF 01761 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01762 extension=".epLog",message=routineP//& 01763 " preconditioner setup main_p_env finished",error=error) 01764 01765 CALL section_vals_val_get(ep_env%root_section,"FORCE_EVAL%EP%E0_ONLY",& 01766 l_val=e0_only,error=error) 01767 IF (.NOT. e0_only) THEN 01768 ! *** perturb to do 01769 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01770 extension=".epLog",message=routineP//& 01771 " perturbation start",error=error) 01772 CALL section_vals_val_get(ep_env%root_section,"FORCE_EVAL%EP%EPS_LIN_SOLV",& 01773 r_val=eps_lin_solv,error=error) 01774 CALL stupid_solve(ep_env,eps_r=eps_lin_solv,error=error) 01775 01776 ep_env%energy%e1=0.0_dp 01777 DO ispin=1,ep_env%nspins 01778 CALL cp_fm_trace(matrix_a=ep_env%m_pi_Hrho_psi0d(ispin)%matrix,& 01779 matrix_b=ep_env%psi1(ispin)%matrix,trace=e_corr_s,error=error) 01780 ep_env%energy%e1=ep_env%energy%e1-e_corr_s 01781 END DO 01782 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01783 extension=".epLog",message=routineP//& 01784 " perturbation end",error=error) 01785 ! *** end perturb 01786 END IF 01787 01788 ep_env%energy%e_tot=ep_env%energy%e0+ep_env%energy%e1 01789 output_unit=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%RUN_INFO",& 01790 extension=".epLog",error=error) 01791 IF (output_unit>0) THEN 01792 WRITE (UNIT=output_unit,FMT="(T3,A,T56,F25.14)")& 01793 "no interaction energy:",ep_env%energy%e_no_int 01794 WRITE (UNIT=output_unit,FMT="(T3,A,T56,F25.14)")& 01795 "psi0 energy (e0):",ep_env%energy%e0 01796 WRITE (UNIT=output_unit,FMT="(T3,A,T56,F25.14)")& 01797 "correction energy (e1):",ep_env%energy%e1 01798 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01799 "Total energy: ",ep_env%energy%e_tot 01800 CALL cp_print_key_finished_output(unit_nr,logger,ep_env%input,& 01801 "PRINT%RUN_INFO",& 01802 error=error) 01803 END IF 01804 01805 output_unit=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%ENERGY",& 01806 extension=".epEnergy",error=error) 01807 IF (output_unit>0) THEN 01808 WRITE (UNIT=output_unit,FMT="(4F25.14)")& 01809 ep_env%energy%e_no_int,ep_env%energy%e0,ep_env%energy%e1,ep_env%energy%e_tot 01810 CALL cp_print_key_finished_output(output_unit,logger,ep_env%input,& 01811 "PRINT%ENERGY", error=error) 01812 END IF 01813 01814 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01815 ep_env%input,"PRINT%EP_MATRIXES/PSI1",error=error),cp_p_file)) THEN 01816 unit_nr=cp_print_key_unit_nr(logger,ep_env%input,& 01817 "PRINT%EP_MATRIXES/PSI1",extension=".psi1",& 01818 error=error) 01819 CALL cp_fm_vect_write(ep_env%psi1, unit_nr,& 01820 long_description=.TRUE., error=suberror) 01821 CALL cp_error_reset(suberror) 01822 CALL cp_print_key_finished_output(unit_nr,logger,ep_env%input,& 01823 "PRINT%EP_MATRIXES/PSI1",& 01824 error=error) 01825 END IF 01826 IF (ASSOCIATED(ep_env%main_p_env%rho1)) THEN 01827 IF (ep_env%main_p_env%rho1%rho_r_valid) THEN 01828 IF (BTEST(cp_print_key_should_output(logger%iter_info,ep_env%input,& 01829 "PRINT%EP_RHO_CUBE/RHO1",error=error),cp_p_file) ) THEN 01830 DO ispin=1,dft_control%nspins 01831 unit_nr=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%EP_RHO_CUBE/RHO1",& 01832 extension=".cube",middle_name="rho1_"//TRIM(ADJUSTL(cp_to_string(ispin))),& 01833 log_filename=.FALSE.,error=error) 01834 CALL pw_to_cube ( ep_env%main_p_env%rho1%rho_r(1)%pw, unit_nr=unit_nr,& 01835 title="rho1_"//TRIM(ADJUSTL(cp_to_string(ispin))),& 01836 particles=particles,& 01837 stride=section_get_ivals(ep_env%input,"PRINT%EP_RHO_CUBE%STRIDE",error=error),& 01838 error=suberror) 01839 CALL cp_error_reset(suberror) 01840 CALL cp_print_key_finished_output(unit_nr,logger,ep_env%input,& 01841 "PRINT%EP_RHO_CUBE/RHO1",error=error) 01842 END DO 01843 END IF 01844 END IF 01845 END IF 01846 01847 IF (BTEST(cp_print_key_should_output(logger%iter_info,ep_env%input,& 01848 "PRINT%EP_RHO_CUBE/RHO0",error=error),cp_p_file) ) THEN 01849 DO ispin=1,dft_control%nspins 01850 unit_nr=cp_print_key_unit_nr(logger,ep_env%input,"PRINT%EP_RHO_CUBE/RHO0",& 01851 extension=".cube",middle_name="rho0_"//TRIM(ADJUSTL(cp_to_string(ispin))),& 01852 log_filename=.FALSE.,error=error) 01853 CALL pw_to_cube ( ep_env%main_qs_env%rho%rho_r(1)%pw, unit_nr=unit_nr,& 01854 title="rho0_"//TRIM(ADJUSTL(cp_to_string(ispin))),& 01855 particles=particles,& 01856 stride=section_get_ivals(ep_env%input,"PRINT%EP_RHO_CUBE%STRIDE",error=error),& 01857 error=suberror) 01858 CALL cp_error_reset(suberror) 01859 CALL cp_print_key_finished_output(unit_nr,logger,ep_env%input,& 01860 "PRINT%EP_RHO_CUBE/RHO0",error=error) 01861 END DO 01862 END IF 01863 01864 CALL get_qs_env(qs_env=qs_env,& 01865 ks_env=ks_env,& 01866 matrix_ks=matrix_ks,& 01867 matrix_s=matrix_s,& 01868 mos=mos,& 01869 rho=rho,error=error) 01870 01871 ! forces 01872 ! Build W matrix (from the orthonormality constraint lambda C^T S C) 01873 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01874 extension=".epLog",message=routineP//& 01875 " starting forces",error=error) 01876 01877 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01878 extension=".epLog",message=routineP//& 01879 " calculate dc/dr forces",error=error) 01880 ! commented out (fawzi 24-Apr-07) 01881 ! CALL ep_calc_f_sub(ep_env,error=error) 01882 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01883 extension=".epLog",message=routineP//& 01884 " finished calculating forces dc/dr",error=error) 01885 01886 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01887 extension=".epLog",message=routineP//& 01888 " creating w_matrix main_qs_env",error=error) 01889 nspin = SIZE(mos) 01890 CALL cp_dbcsr_allocate_matrix_set(matrix_w,nspin,error=error) 01891 DO ispin=1,nspin 01892 ALLOCATE(matrix_w(ispin)%matrix) 01893 CALL cp_dbcsr_init(matrix_w(ispin)%matrix,error=error) 01894 CALL cp_dbcsr_copy(matrix_w(ispin)%matrix,matrix_s(1)%matrix,& 01895 "W MATRIX"//TRIM(ADJUSTL(cp_to_string(ispin))),error=error) 01896 END DO 01897 01898 CALL ep_calc_w_matrix_full(ep_env,matrix_w,error=error) 01899 01900 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01901 qs_env%input,"DFT%PRINT%AO_MATRICES/W_MATRIX",error=error),cp_p_file)) THEN 01902 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/W_MATRIX",& 01903 extension=".Log",error=error) 01904 !FM DO ispin=1,nspin 01905 !FM CALL write_sparse_matrix(matrix_w(ispin)%matrix,4,6,qs_env,globenv,output_unit=iw) 01906 !FM END DO 01907 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 01908 "DFT%PRINT%AO_MATRICES/W_MATRIX", error=error) 01909 END IF 01910 01911 CALL set_qs_env(qs_env=qs_env,matrix_w=matrix_w,error=error) 01912 01913 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01914 extension=".epLog",message=routineP//& 01915 " w_matrix created and added",error=error) 01916 01917 !*** compute core forces (also overwrites matrix_w) *** 01918 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01919 extension=".epLog",message=routineP//& 01920 " compute forces in main_qs_env",error=error) 01921 CALL build_core_hamiltonian_matrix(qs_env=qs_env,calculate_forces=.TRUE.,error=error) 01922 CALL calculate_ecore_self(qs_env,error=error) 01923 CALL calculate_ecore_overlap(qs_env, para_env, & 01924 calculate_forces=.TRUE.,& 01925 error=error) 01926 01927 ! *** compute grid-based forces *** 01928 CALL qs_ks_update_qs_env(ks_env=ks_env,& 01929 qs_env=qs_env,& 01930 calculate_forces=.TRUE.,error=error) 01931 01932 ! *** replicate forces *** 01933 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01934 extension=".epLog",message=routineP//& 01935 " summing up forces",error=error) 01936 DO ikind=1,SIZE(force) 01937 CALL mp_sum(force(ikind)%overlap,para_env%group) 01938 CALL mp_sum(force(ikind)%kinetic,para_env%group) 01939 CALL mp_sum(force(ikind)%gth_ppl,para_env%group) 01940 CALL mp_sum(force(ikind)%gth_nlcc,para_env%group) 01941 CALL mp_sum(force(ikind)%gth_ppnl,para_env%group) 01942 CALL mp_sum(force(ikind)%all_potential,para_env%group) 01943 CALL mp_sum(force(ikind)%core_overlap,para_env%group) 01944 CALL mp_sum(force(ikind)%rho_core,para_env%group) 01945 CALL mp_sum(force(ikind)%rho_elec,para_env%group) 01946 CALL mp_sum(force(ikind)%vhxc_atom,para_env%group) 01947 CALL mp_sum(force(ikind)%g0s_Vh_elec,para_env%group) 01948 CALL mp_sum(force(ikind)%fock_4c,para_env%group) 01949 force(ikind)%total(:,:) = force(ikind)%total(:,:) +& 01950 force(ikind)%core_overlap(:,:) +& 01951 force(ikind)%gth_ppl(:,:) +& 01952 force(ikind)%gth_nlcc(:,:) +& 01953 force(ikind)%gth_ppnl(:,:) +& 01954 force(ikind)%all_potential(:,:) +& 01955 force(ikind)%kinetic(:,:) +& 01956 force(ikind)%overlap(:,:) +& 01957 force(ikind)%rho_core(:,:) +& 01958 force(ikind)%rho_elec(:,:) +& 01959 force(ikind)%vhxc_atom(:,:) +& 01960 force(ikind)%g0s_Vh_elec(:,:) +& 01961 force(ikind)%fock_4c(:,:) 01962 END DO 01963 01964 DO iatom=1,natom 01965 01966 i = atom_of_kind(iatom) 01967 ikind = kind_of(iatom) 01968 ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 01969 ! the force is - dE/dR, what is called force is actually the gradient 01970 ! Things should have the right name 01971 ! The minus sign below is a hack 01972 ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 01973 force(ikind)%other(1:3,i)=-particle_set(iatom)%f(1:3)+force(ikind)%ch_pulay(1:3,i)+ep_env%force%f0_internal(1:3,iatom) 01974 force(ikind)%total(1:3,i)=force(ikind)%total(1:3,i)+force(ikind)%other(1:3,i) 01975 particle_set(iatom)%f = -force(ikind)%total(1:3,i) 01976 END DO 01977 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 01978 extension=".epLog",message=routineP//& 01979 " finished forces",error=error) 01980 01981 output_unit = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%DERIVATIVES",& 01982 extension=".Log",error=error) 01983 print_section => section_vals_get_subs_vals(qs_env%input,"DFT%PRINT%DERIVATIVES",error=error) 01984 CALL write_forces(force,atomic_kind_set,0,output_unit=output_unit,& 01985 print_section=print_section,error=error) 01986 01987 CALL cp_print_key_finished_output(output_unit,logger,qs_env%input,& 01988 "DFT%PRINT%DERIVATIVES",error=error) 01989 01990 ! destroy preconditioner 01991 IF (ASSOCIATED(ep_env%precond)) THEN 01992 !CASE(ot_precond_full_all,ot_precond_full_single, ot_precond_full_single_inverse) 01993 ! these depend on the ks matrix 01994 DO ispin=1,SIZE(ep_env%precond) 01995 CALL destroy_preconditioner(ep_env%precond(ispin)%preconditioner,error=error) 01996 DEALLOCATE(ep_env%precond(ispin)%preconditioner) 01997 ENDDO 01998 DEALLOCATE(ep_env%precond,stat=stat) 01999 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02000 !CASE(ot_precond_none,ot_precond_full_kinetic,ot_precond_s_inverse, & 02001 ! ot_precond_sparse_diag) ! these are 'independent' 02002 END IF 02003 02004 CALL cp_dbcsr_deallocate_matrix_set( matrix_w, error=error ) 02005 02006 CALL set_qs_env(qs_env=qs_env,matrix_w=matrix_w,error=error) 02007 02008 DEALLOCATE (atom_of_kind,STAT=stat) 02009 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 02010 02011 DEALLOCATE (kind_of,STAT=stat) 02012 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 02013 02014 CALL cp_error_dealloc_ref(suberror) 02015 CALL cp_print_key_log(logger,ep_env%input,"PRINT%RUN_INFO",& 02016 extension=".epLog",message=routineP//& 02017 " finished",error=error) 02018 02019 CALL timestop(handle) 02020 END SUBROUTINE ep_env_calc_e0 02021 02022 ! ***************************************************************************** 02028 SUBROUTINE ep_calc_w_matrix_full(ep_env,w_matrix,error) 02029 TYPE(ep_env_type), POINTER :: ep_env 02030 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02031 POINTER :: w_matrix 02032 TYPE(cp_error_type), INTENT(inout) :: error 02033 02034 CHARACTER(len=*), PARAMETER :: routineN = 'ep_calc_w_matrix_full', 02035 routineP = moduleN//':'//routineN 02036 02037 INTEGER :: handle, ierr, ispin 02038 LOGICAL :: failure 02039 TYPE(cp_error_type) :: new_error 02040 TYPE(cp_fm_pool_type), POINTER :: maxao_maxmo_fm_pool 02041 TYPE(cp_fm_type), POINTER :: psi0, weighted_vectors 02042 TYPE(f_env_type), POINTER :: f_env 02043 TYPE(mo_set_p_type), DIMENSION(:), 02044 POINTER :: mos 02045 02046 CALL timeset(routineN,handle) 02047 failure=.FALSE. 02048 NULLIFY(weighted_vectors,f_env,maxao_maxmo_fm_pool,mos) 02049 02050 CALL f_env_add_defaults(ep_env%f_env_id,f_env,new_error,failure) 02051 CALL mpools_get(ep_env%main_qs_env%mpools,& 02052 maxao_maxmo_fm_pool=maxao_maxmo_fm_pool,& 02053 error=error) 02054 CALL fm_pool_create_fm(maxao_maxmo_fm_pool,weighted_vectors,"weighted_vectors",& 02055 error=error) 02056 CALL get_qs_env(ep_env%main_qs_env,mos=mos,error=new_error) 02057 DO ispin=1,ep_env%nspins 02058 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=psi0) 02059 CALL cp_fm_symm('R', 'U', ep_env%full_nao(ispin), ep_env%full_nmo(ispin), & 02060 -REAL(3-ep_env%nspins,dp), 02061 ep_env%main_p_env%m_epsilon(ispin)%matrix,psi0, 0.0_dp,weighted_vectors,error=error) 02062 CALL cp_dbcsr_set(w_matrix(ispin)%matrix,0.0_dp,error=error) 02063 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix(ispin)%matrix,& 02064 matrix_v=psi0,& 02065 matrix_g=weighted_vectors,& 02066 ncol=ep_env%full_nmo(ispin),error=error) 02067 END DO 02068 CALL fm_pool_give_back_fm(maxao_maxmo_fm_pool,weighted_vectors,error=error) 02069 CALL f_env_rm_defaults(f_env,new_error,ierr) 02070 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 02071 CALL timestop(handle) 02072 02073 END SUBROUTINE ep_calc_w_matrix_full 02074 02075 ! ***************************************************************************** 02089 SUBROUTINE ep_calc_dc_dr_fdiff(ep_env,dc_dr,imol,iat,idir,error) 02090 TYPE(ep_env_type), POINTER :: ep_env 02091 TYPE(cp_fm_p_type), DIMENSION(:) :: dc_dr 02092 INTEGER, INTENT(in) :: imol, iat, idir 02093 TYPE(cp_error_type), INTENT(inout) :: error 02094 02095 CHARACTER(len=*), PARAMETER :: routineN = 'ep_calc_dc_dr_fdiff', 02096 routineP = moduleN//':'//routineN 02097 REAL(kind=dp), PARAMETER :: dr = 0.001 02098 02099 INTEGER :: handle, ierr, ilocal, ipos, 02100 ispin, ndims 02101 LOGICAL :: failure 02102 REAL(kind=dp) :: e_pot, orig 02103 TYPE(cp_fm_pool_p_type), DIMENSION(:), 02104 POINTER :: mo_mo_fm_pools 02105 TYPE(cp_fm_p_type), DIMENSION(:), 02106 POINTER :: mo_mo_m 02107 TYPE(cp_error_type) :: new_error 02108 REAL(kind=dp), 02109 DIMENSION(ep_env%mol_envs%ndim) :: pos 02110 TYPE(f_env_type), POINTER :: f_env 02111 TYPE(qs_environment_type), POINTER :: sub_qs_env 02112 TYPE(qs_matrix_pools_type), POINTER :: mpools 02113 TYPE(qs_wf_history_type), POINTER :: tmp_wf_history 02114 TYPE(qs_wf_snapshot_type), POINTER :: snapshot, tmp_snap 02115 TYPE(replica_env_type), POINTER :: rep_env 02116 02117 failure=.FALSE. 02118 02119 CALL timeset(routineN,handle) 02120 NULLIFY(snapshot,tmp_snap,mo_mo_m) 02121 CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 02122 CPPrecondition(ep_env%ref_count>0,cp_failure_level,routineP,error,failure) 02123 IF (.NOT.failure) THEN 02124 rep_env => ep_env%mol_envs 02125 CALL f_env_add_defaults(rep_env%f_env_id,f_env,new_error,failure) 02126 CALL force_env_get(f_env%force_env,qs_env=sub_qs_env,& 02127 error=new_error) 02128 CALL get_qs_env(sub_qs_env,mpools=mpools,error=new_error) 02129 CALL mpools_get(mpools,mo_mo_fm_pools=mo_mo_fm_pools,error=new_error) 02130 02131 ! ** create tmp wf history 02132 ilocal=rep_env_local_index(rep_env,imol,error=error) 02133 snapshot => wfi_get_snapshot(& 02134 rep_env%wf_history(ilocal)%wf_history,& 02135 1,error=new_error) 02136 CPAssert(ASSOCIATED(snapshot),cp_failure_level,routineP,new_error,failure) 02137 CPAssert(ASSOCIATED(snapshot%wf),cp_failure_level,routineP,new_error,failure) 02138 CALL wfi_create(tmp_wf_history,& 02139 interpolation_method_nr=wfi_use_prev_wf_method_nr,& 02140 extrapolation_order = 1,& 02141 has_unit_metric = .FALSE.,& 02142 error=new_error) 02143 CALL wfi_change_memory_depth(tmp_wf_history,& 02144 MAX(1,tmp_wf_history%memory_depth),error=error) 02145 tmp_wf_history%store_wf=.TRUE. 02146 02147 CALL wfs_duplicate_snapshot(input_snapshot=snapshot,& 02148 output_snapshot=tmp_snap,& 02149 qs_env=sub_qs_env, error=new_error) 02150 tmp_wf_history%last_state_index=MODULO(tmp_wf_history%snapshot_count,& 02151 tmp_wf_history%memory_depth)+1 02152 CALL wfs_release(tmp_wf_history%past_states(tmp_wf_history%last_state_index)%snapshot,& 02153 error=new_error) 02154 tmp_wf_history%past_states(tmp_wf_history%last_state_index)%snapshot => & 02155 tmp_snap 02156 NULLIFY(tmp_snap) 02157 CALL set_qs_env(sub_qs_env,wf_history=tmp_wf_history,error=new_error) 02158 02159 ndims=ep_env%mol_envs%ndim 02160 pos=rep_env%r(:,imol) 02161 IF (.not.failure) THEN 02162 ipos=3*(iat-1)+idir 02163 orig=pos(ipos) 02164 pos(ipos)=orig+dr 02165 CALL calc_energy(env_id=ep_env%mol_envs%f_env_id,pos=pos,n_el=ndims,& 02166 e_pot=e_pot,ierr=ierr) 02167 02168 END IF 02169 IF (.NOT.failure) THEN 02170 tmp_snap => wfi_get_snapshot(& 02171 tmp_wf_history,& 02172 1,error=new_error) 02173 CPAssert(ASSOCIATED(tmp_snap),cp_failure_level,routineP,new_error,failure) 02174 CPAssert(ASSOCIATED(tmp_snap%wf),cp_failure_level,routineP,new_error,failure) 02175 CALL fm_pools_create_fm_vect(mo_mo_fm_pools,mo_mo_m,error=new_error) 02176 DO ispin=1,ep_env%nspins 02177 CALL cp_dbcsr_sm_fm_multiply(snapshot%overlap,& 02178 tmp_snap%wf(ispin)%matrix,dc_dr(ispin)%matrix,& 02179 ncol=ep_env%sub_nmo(ispin),error=new_error) 02180 CALL cp_fm_gemm(transa='T',transb='N',m=ep_env%sub_nmo(ispin),& 02181 n=ep_env%sub_nmo(ispin),k=ep_env%sub_nao(ispin),alpha=1.0_dp,& 02182 matrix_a=snapshot%wf(ispin)%matrix,& 02183 matrix_b=dc_dr(ispin)%matrix,beta=0._dp,& 02184 matrix_c=mo_mo_m(ispin)%matrix,error=new_error) 02185 CALL cp_fm_scale_and_add(matrix_a=dc_dr(ispin)%matrix,& 02186 matrix_b=tmp_snap%wf(ispin)%matrix,alpha=0.0_dp,beta=1.0_dp,& 02187 error=new_error) 02188 CALL cp_fm_gemm(transa='N',transb='N',m=ep_env%sub_nao(ispin),& 02189 n=ep_env%sub_nmo(ispin),k=ep_env%sub_nmo(ispin),alpha=-1.0_dp/dr,& 02190 matrix_a=snapshot%wf(ispin)%matrix,& 02191 matrix_b=mo_mo_m(ispin)%matrix,beta=1._dp/dr,& 02192 matrix_c=dc_dr(ispin)%matrix,error=new_error) 02193 END DO 02194 CALL fm_pools_give_back_fm_vect(mo_mo_fm_pools,mo_mo_m,error=new_error) 02195 END IF 02196 CALL wfi_release(tmp_wf_history,error=new_error) 02197 02198 CALL f_env_rm_defaults(f_env,new_error,ierr) 02199 CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 02200 END IF 02201 CALL timestop(handle) 02202 END SUBROUTINE ep_calc_dc_dr_fdiff 02203 02204 ! ***************************************************************************** 02205 ! The following routines are never called, but introduce a dependence on the 02206 ! old SM code. 02214 02215 ! CHARACTER(len=*), PARAMETER :: routineN = 'ep_calc_f_sub', & 02216 ! routineP = moduleN//':'//routineN 02217 02218 ! INTEGER :: handle, iat, idir, ierr, & 02219 ! irep, irep_local, ispin, nat, & 02220 ! stat 02221 ! INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 02222 ! LOGICAL :: failure 02223 ! REAL(kind=dp) :: fdiff(2) 02224 ! TYPE(atomic_kind_type), DIMENSION(:), & 02225 ! POINTER :: atomic_kind_set 02226 ! TYPE(cp_dbcsr_p_type), DIMENSION(:), & 02227 ! POINTER :: matrix_s_b 02228 ! TYPE(cp_error_type) :: new_error 02229 ! TYPE(cp_fm_p_type), DIMENSION(:), & 02230 ! POINTER :: aomo_work1, coeffs_cr, & 02231 ! pi_Hrho_psi0d 02232 ! TYPE(cp_fm_pool_p_type), DIMENSION(:), & 02233 ! POINTER :: ao_mo_fm_pools 02234 ! TYPE(cp_fm_pool_type), POINTER :: maxao_maxao_fm_pool, & 02235 ! maxao_maxmo_fm_pool, & 02236 ! maxmo_maxmo_fm_pool 02237 ! TYPE(cp_fm_type), POINTER :: aoao_work1, psi0 02238 ! TYPE(ep_force_type), POINTER :: ep_force 02239 ! TYPE(f_env_type), POINTER :: f_env 02240 ! TYPE(global_environment_type), POINTER :: sub_globenv 02241 ! TYPE(mo_set_p_type), DIMENSION(:), & 02242 ! POINTER :: mos 02243 ! TYPE(qs_environment_type), POINTER :: sub_qs_env 02244 ! TYPE(qs_force_type), DIMENSION(:), & 02245 ! POINTER :: force 02246 ! TYPE(qs_wf_history_type), POINTER :: tmp_wf_history 02247 ! TYPE(qs_wf_snapshot_type), POINTER :: snapshot, tmp_snap 02248 ! TYPE(real_matrix_p_type), DIMENSION(:), & 02249 ! POINTER :: matrix_s 02250 ! TYPE(real_matrix_p_type), & 02251 ! DIMENSION(:, :), POINTER :: local_P_H_pi 02252 ! TYPE(real_matrix_type), POINTER :: diag_P_H_pi 02253 ! TYPE(replica_env_type), POINTER :: rep_env 02254 02255 ! CALL timeset(routineN,handle) 02256 ! failure=.FALSE. 02257 ! NULLIFY(f_env, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool, & 02258 ! maxao_maxmo_fm_pool, ao_mo_fm_pools, mos, sub_globenv, sub_qs_env, distribution_2d) 02259 ! ep_force => ep_env%force 02260 02261 ! CALL get_qs_env(ep_env%main_qs_env,mos=mos,matrix_s=matrix_s_b,& 02262 ! distribution_2d=distribution_2d,& 02263 ! error=new_error) 02264 02265 ! NULLIFY(matrix_s)!sm->dbcsr 02266 ! CALL allocate_matrix_set( matrix_s, SIZE(matrix_s_b), error )!sm->dbcsr 02267 ! DO ispin=1,SIZE(matrix_s)!sm->dbcsr 02268 ! CALL sm_from_dbcsr(matrix_s(ispin)%matrix, matrix_s_b(ispin)%matrix, & 02269 ! distribution_2d,error)!sm->dbcsr 02270 ! ENDDO!sm->dbcsr 02271 02272 ! NULLIFY(diag_P_H_pi,local_P_H_pi,local_P_H_pi) 02273 ! CALL ep_create_local_empty_sm(ep_env,local_P_H_pi,"local_P_H_pi",& 02274 ! symmetry="none",error=error) 02275 ! CALL ep_create_full_diag_sm(ep_env,matrix_s(1)%matrix,diag_P_H_pi,"diag_P_H_pi",error=error) 02276 ! DO ispin=1,ep_env%nspins 02277 ! CALL set_matrix(diag_P_H_pi,0.0_dp) 02278 ! CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=psi0) 02279 ! CALL cp_sm_plus_fm_fm_t(sparse_matrix=diag_P_H_pi,& 02280 ! matrix_v=psi0,& 02281 ! matrix_g=ep_env%m_pi_Hrho_psi0d(ispin)%matrix,& 02282 ! ncol=ep_env%full_nmo(ispin),error=error) 02283 !FM PRINT *,"diag_P_H_pi" 02284 !FM CALL write_sparse_matrix(diag_P_H_pi,4,6,ep_env%main_qs_env,ep_env%globenv,output_unit=6) 02285 ! CALL ep_full_to_local_sm(ep_env,diag_P_H_pi,local_P_H_pi,ispin,& 02286 ! error=error) 02287 ! END DO 02288 ! CALL deallocate_matrix(diag_P_H_pi,error=error) 02289 02290 ! rep_env => ep_env%mol_envs 02291 ! CALL f_env_add_defaults(rep_env%f_env_id,f_env,new_error,failure) 02292 ! CALL force_env_get(f_env%force_env,qs_env=sub_qs_env,globenv=sub_globenv,& 02293 ! error=new_error) 02294 !FM DO ispin=1,ep_env%nspins 02295 !FM DO i=1,SIZE(local_P_H_pi,2) 02296 !FM PRINT *,"local_P_H_pi",i 02297 !FM CALL write_sparse_matrix(local_P_H_pi(ispin,i)%matrix,& 02298 !FM 4,6,sub_qs_env,sub_globenv,output_unit=6) 02299 !FM END DO 02300 !FM END DO 02301 ! CALL mpools_get(sub_qs_env%mpools,& 02302 ! maxao_maxao_fm_pool=maxao_maxao_fm_pool,& 02303 ! maxao_maxmo_fm_pool=maxao_maxmo_fm_pool,& 02304 ! maxmo_maxmo_fm_pool=maxao_maxmo_fm_pool,& 02305 ! ao_mo_fm_pools=ao_mo_fm_pools,& 02306 ! error=error) 02307 ! CALL fm_pools_create_fm_vect(ao_mo_fm_pools,aomo_work1,& 02308 ! name="ao_mo_work1",error=error) 02309 ! CALL fm_pools_create_fm_vect(ao_mo_fm_pools,coeffs_cr,& 02310 ! name="coeffs_cr",error=error) 02311 02312 ! CALL get_qs_env(sub_qs_env,force=force,atomic_kind_set=atomic_kind_set,& 02313 ! error=error) 02314 ! nat=rep_env%nat 02315 ! ALLOCATE (atom_of_kind(nat),STAT=stat) 02316 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02317 ! ALLOCATE (kind_of(nat),STAT=stat) 02318 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02319 ! CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 02320 ! atom_of_kind=atom_of_kind,& 02321 ! kind_of=kind_of) 02322 02323 ! DO irep_local=1,SIZE(rep_env%local_rep_indices) 02324 ! irep=rep_env%local_rep_indices(irep_local) 02325 ! snapshot => wfi_get_snapshot(& 02326 ! rep_env%wf_history(irep_local)%wf_history,& 02327 ! 1,error=new_error) 02328 02329 ! CPAssert(ASSOCIATED(snapshot),cp_failure_level,routineP,new_error,failure) 02330 ! CPAssert(ASSOCIATED(snapshot%wf),cp_failure_level,routineP,new_error,failure) 02331 ! CALL get_qs_env(sub_qs_env,error=new_error) 02332 02333 ! NULLIFY(tmp_wf_history,tmp_snap) 02334 ! CALL wfi_create(tmp_wf_history,& 02335 ! interpolation_method_nr=wfi_use_prev_wf_method_nr,& 02336 ! extrapolation_order = 1,& 02337 ! error=new_error) 02338 ! CALL wfi_change_memory_depth(tmp_wf_history,MAX(1,tmp_wf_history%memory_depth),& 02339 ! error=error) 02340 ! tmp_wf_history%store_wf=.TRUE. 02341 ! CALL wfs_duplicate_snapshot(input_snapshot=snapshot, output_snapshot=tmp_snap,& 02342 ! qs_env=sub_qs_env, error=new_error) 02343 ! tmp_wf_history%last_state_index=MODULO(tmp_wf_history%snapshot_count,& 02344 ! tmp_wf_history%memory_depth)+1 02345 ! CALL wfs_release(tmp_wf_history%past_states(tmp_wf_history%last_state_index)%snapshot,& 02346 ! error=new_error) 02347 ! tmp_wf_history%past_states(tmp_wf_history%last_state_index)%snapshot => & 02348 ! tmp_snap 02349 ! NULLIFY(tmp_snap) 02350 ! CALL set_qs_env(sub_qs_env,wf_history=tmp_wf_history,error=new_error) 02351 ! CALL wfi_release(tmp_wf_history,error=error) 02352 02353 ! !* calc H_0 subforces and deriv operators 02354 ! CALL calc_energy_force(rep_env%f_env_id,.TRUE.,ierr) 02355 ! CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 02356 02357 ! !* setup sub perturb env 02358 ! NULLIFY(pi_Hrho_psi0d) 02359 ! CALL fm_pools_create_fm_vect(ao_mo_fm_pools,pi_Hrho_psi0d,error=error) 02360 ! CALL p_env_did_change(ep_env%sub_p_env, s_struct_changed=.TRUE., & 02361 ! grid_changed=.FALSE.,error=new_error) 02362 ! CALL p_env_psi0_changed(ep_env%sub_p_env, qs_env=sub_qs_env,& 02363 ! Hrho_psi0d=pi_Hrho_psi0d, error=new_error) 02364 ! CALL p_postortho(p_env=ep_env%sub_p_env,qs_env=sub_qs_env,& 02365 ! v=pi_Hrho_psi0d, & 02366 ! n_cols=ep_env%sub_p_env%n_mo,& 02367 ! error=error) 02368 02369 ! !* calc dC/dR coefficients 02370 ! CALL fm_pool_create_fm(maxao_maxao_fm_pool,aoao_work1,& 02371 ! name="ao_ao_work1",error=error) 02372 ! DO ispin=1,ep_env%nspins 02373 ! CALL copy_sm_to_fm(local_P_H_pi(ispin,irep_local)%matrix,aoao_work1,error=new_error) 02374 ! CALL cp_fm_gemm(transa="T",transb="N",m=ep_env%sub_nao(ispin),& 02375 ! n=ep_env%sub_nmo(ispin),k=ep_env%sub_nao(ispin),alpha=1._dp,& 02376 ! matrix_a=aoao_work1,matrix_b=ep_env%sub_p_env%S_psi0(ispin)%matrix,& 02377 ! beta=0._dp,matrix_c=coeffs_cr(ispin)%matrix,error=new_error) 02378 ! END DO 02379 ! CALL fm_pool_give_back_fm(maxao_maxao_fm_pool,aoao_work1,error=error) 02380 02381 ! ! invert second order matrix and finds coeffs of ddE/(dC dR) 02382 ! ! analytic 02383 ! !CALL ep_sub_stupid_solve(ep_env,coeffs_cr,coeffs_ecr) 02384 ! ! f_diff 02385 ! DO iat=1,ep_env%nat_per_mol 02386 ! DO idir=1,3 02387 ! ! f_diff deriv 02388 ! CALL ep_calc_dc_dr_fdiff(ep_env,aomo_work1,irep,iat,idir,error) 02389 ! DO ispin=1,ep_env%nspins 02390 ! CALL cp_fm_trace(coeffs_cr(ispin)%matrix,& 02391 ! aomo_work1(ispin)%matrix,fdiff(ispin), error=error) 02392 ! END DO 02393 02394 ! DO ispin=1,ep_env%nspins 02395 ! ep_force%f0_internal(idir,iat)=ep_force%f0_internal(idir,iat)+fdiff(ispin) 02396 !FM PRINT *,"irep",irep," iatom",iatom 02397 !FM PRINT *," ep_force%f0_internal(1:3,atoms(iatom))",ep_force%f0_internal(1:3,atoms(iatom)) 02398 ! END DO 02399 ! END DO 02400 ! END DO 02401 ! CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,pi_Hrho_psi0d,error=error) 02402 ! END DO 02403 ! CALL ep_give_back_local_sm(ep_env,local_P_H_pi,error=error) 02404 ! CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,aomo_work1,error=error) 02405 ! CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools,coeffs_cr,error=error) 02406 02407 ! DEALLOCATE (atom_of_kind,STAT=stat) 02408 ! CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 02409 02410 ! DEALLOCATE (kind_of,STAT=stat) 02411 ! CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 02412 02413 ! CALL f_env_rm_defaults(f_env,new_error,ierr) 02414 ! CPAssert(ierr==0,cp_failure_level,routineP,error,failure) 02415 02416 ! CALL mp_sum(ep_force%f0_internal,ep_env%main_qs_env%para_env%group) 02417 02418 ! CALL deallocate_matrix_set( matrix_s, error )!sm->dbcsr 02419 02420 ! CALL timestop(handle) 02421 02422 ! END SUBROUTINE ep_calc_f_sub 02423 02424 ! ***************************************************************************** 02437 02438 ! ! History: - Creation (17.11.2000, MK) 02439 02440 ! TYPE(ep_env_type), POINTER :: ep_env 02441 ! TYPE(real_matrix_type), POINTER :: source_m, diag_m 02442 ! CHARACTER(LEN=*), INTENT(IN) :: diag_m_name 02443 ! TYPE(cp_error_type), INTENT(inout) :: error 02444 02445 ! CHARACTER(LEN=*), PARAMETER :: routineN = 'ep_create_full_diag_sm', & 02446 ! routineP = moduleN//':'//routineN 02447 02448 ! INTEGER :: iblock_col, iblock_row, & 02449 ! ncols, nrows 02450 ! INTEGER, DIMENSION(:), POINTER :: cols, rows 02451 ! LOGICAL :: failure 02452 02453 ! failure=.FALSE. 02454 ! CPPrecondition(ASSOCIATED(source_m),cp_failure_level,routineP,error,failure) 02455 ! CPPrecondition(.NOT.ASSOCIATED(diag_m),cp_failure_level,routineP,error,failure) 02456 02457 ! PRINT *,"start ", routineP 02458 ! NULLIFY(rows,cols) 02459 ! CALL allocate_matrix(matrix=diag_m,& 02460 ! nrow=source_m%nrow,& 02461 ! ncol=source_m%ncol,& 02462 ! nblock_row=source_m%nblock_row,& 02463 ! nblock_col=source_m%nblock_col,& 02464 ! first_row=source_m%first_row(:),& 02465 ! last_row=source_m%last_row(:),& 02466 ! first_col=source_m%first_col(:),& 02467 ! last_col=source_m%last_col(:),& 02468 ! matrix_name=diag_m_name,& 02469 ! sparsity_id=source_m%sparsity_id, & 02470 ! distribution_2d=source_m%distribution_2d,& 02471 ! matrix_symmetry="none",error=error) 02472 02473 ! CALL distribution_2d_get(source_m%distribution_2d,& 02474 ! flat_local_rows=rows,flat_local_cols=cols,& 02475 ! n_flat_local_rows=nrows,n_flat_local_cols=ncols,& 02476 ! error=error) 02477 02478 ! DO iblock_row=1,nrows 02479 ! DO iblock_col=ncols,1,-1 02480 ! IF (ep_env%at2sub(cols(iblock_col))==ep_env%at2sub(rows(iblock_row))) THEN 02481 ! PRINT *, "adding ",iblock_row,iblock_col 02482 ! CALL add_block_node(matrix=diag_m,& 02483 ! block_row=rows(iblock_row),& 02484 ! block_col=cols(iblock_col),error=error) 02485 ! END IF 02486 ! END DO 02487 ! END DO 02488 ! PRINT *,"end ", routineP 02489 02490 ! END SUBROUTINE ep_create_full_diag_sm 02491 02492 02493 ! ***************************************************************************** 02510 02511 ! CHARACTER(len=*), PARAMETER :: routineN = 'ep_create_local_empty_sm', & 02512 ! routineP = moduleN//':'//routineN 02513 02514 ! INTEGER :: handle, i, ierr, ispin, stat 02515 ! LOGICAL :: failure 02516 ! TYPE(cp_dbcsr_p_type), DIMENSION(:), & 02517 ! POINTER :: matrix_s_b 02518 ! TYPE(cp_error_type) :: new_error 02519 ! TYPE(distribution_2d_type), POINTER :: distribution_2d 02520 ! TYPE(f_env_type), POINTER :: f_env 02521 ! TYPE(qs_environment_type), POINTER :: sub_qs_env 02522 ! TYPE(real_matrix_p_type), DIMENSION(:), & 02523 ! POINTER :: matrix_s 02524 ! TYPE(real_matrix_type), POINTER :: source_m 02525 02526 ! failure=.FALSE. 02527 ! CALL timeset(routineN,handle) 02528 ! NULLIFY(f_env,sub_qs_env,matrix_s,distribution_2d) 02529 ! CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 02530 ! CPPrecondition(ep_env%ref_count>0,cp_failure_level,routineP,error,failure) 02531 ! IF (.NOT. failure) THEN 02532 ! ALLOCATE(local_sm(ep_env%nspins,SIZE(ep_env%mol_envs%local_rep_indices)),& 02533 ! stat=stat) 02534 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02535 ! END IF 02536 ! IF (.not.failure) THEN 02537 ! CALL f_env_add_defaults(ep_env%f_env_id,f_env=f_env,& 02538 ! new_error=new_error,failure=failure) 02539 ! CALL force_env_get(f_env%force_env,qs_env=sub_qs_env,error=new_error) 02540 02541 ! CALL get_qs_env(sub_qs_env,matrix_s=matrix_s_b,& 02542 ! distribution_2d=distribution_2d,& 02543 ! error=new_error) 02544 02545 ! NULLIFY(matrix_s)!sm->dbcsr 02546 ! CALL allocate_matrix_set( matrix_s, SIZE(matrix_s_b), error )!sm->dbcsr 02547 ! DO ispin=1,SIZE(matrix_s)!sm->dbcsr 02548 ! CALL sm_from_dbcsr(matrix_s(ispin)%matrix, matrix_s_b(ispin)%matrix, & 02549 ! distribution_2d,error)!sm->dbcsr 02550 ! ENDDO!sm->dbcsr 02551 02552 02553 ! DO i=1,SIZE(local_sm,2) 02554 ! DO ispin=1,ep_env%nspins 02555 ! source_m => matrix_s(ispin)%matrix 02556 ! NULLIFY(local_sm(ispin,i)%matrix) 02557 ! IF (PRESENT(symmetry)) THEN 02558 ! CALL allocate_matrix(matrix=local_sm(ispin,i)%matrix,& 02559 ! nrow=source_m%nrow,& 02560 ! ncol=source_m%ncol,& 02561 ! nblock_row=source_m%nblock_row,& 02562 ! nblock_col=source_m%nblock_col,& 02563 ! first_row=source_m%first_row,& 02564 ! last_row=source_m%last_row,& 02565 ! first_col=source_m%first_col,& 02566 ! last_col=source_m%last_col,& 02567 ! matrix_name=name//TRIM(ADJUSTL(cp_to_string(ispin)))//"_"& 02568 ! //TRIM(ADJUSTL(cp_to_string(i))),& 02569 ! sparsity_id=source_m%sparsity_id+100, & 02570 ! distribution_2d=source_m%distribution_2d,& 02571 ! matrix_symmetry=symmetry,error=error) 02572 ! ELSE 02573 ! CALL allocate_matrix(matrix=local_sm(ispin,i)%matrix,& 02574 ! nrow=source_m%nrow,& 02575 ! ncol=source_m%ncol,& 02576 ! nblock_row=source_m%nblock_row,& 02577 ! nblock_col=source_m%nblock_col,& 02578 ! first_row=source_m%first_row,& 02579 ! last_row=source_m%last_row,& 02580 ! first_col=source_m%first_col,& 02581 ! last_col=source_m%last_col,& 02582 ! matrix_name=name//TRIM(ADJUSTL(cp_to_string(ispin)))//"_"& 02583 ! //TRIM(ADJUSTL(cp_to_string(i))),& 02584 ! sparsity_id=source_m%sparsity_id+100, & 02585 ! distribution_2d=source_m%distribution_2d,& 02586 ! matrix_symmetry=source_m%symmetry,error=error) 02587 ! END IF 02588 ! END DO 02589 ! END DO 02590 ! CALL f_env_rm_defaults(f_env,new_error,ierr) 02591 ! CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 02592 ! END IF 02593 02594 ! DO ispin=1,SIZE(matrix_s)!sm->dbcsr 02595 ! CALL cp_dbcsr_from_sm(matrix_s_b(ispin)%matrix, matrix_s(ispin)%matrix, error)!sm->dbcsr 02596 ! ENDDO!sm->dbcsr 02597 ! CALL deallocate_matrix_set( matrix_s, error )!sm->dbcsr 02598 02599 ! CALL timestop(handle) 02600 ! END SUBROUTINE ep_create_local_empty_sm 02601 02602 ! ***************************************************************************** 02615 02616 ! CHARACTER(len=*), PARAMETER :: routineN = 'ep_give_back_local_sm', & 02617 ! routineP = moduleN//':'//routineN 02618 02619 ! INTEGER :: i, ierr, ispin, stat 02620 ! LOGICAL :: failure 02621 ! TYPE(cp_error_type) :: new_error 02622 ! TYPE(f_env_type), POINTER :: f_env 02623 ! TYPE(qs_environment_type), POINTER :: sub_qs_env 02624 02625 ! failure=.FALSE. 02626 ! NULLIFY(f_env,sub_qs_env) 02627 ! CPPrecondition(ASSOCIATED(ep_env),cp_failure_level,routineP,error,failure) 02628 ! CPPrecondition(ep_env%ref_count>0,cp_failure_level,routineP,error,failure) 02629 ! CPPrecondition(ASSOCIATED(local_sm),cp_failure_level,routineP,error,failure) 02630 ! CPPrecondition(SIZE(local_sm,1)==ep_env%nspins,cp_failure_level,routineP,error,failure) 02631 ! IF (.NOT. failure) THEN 02632 ! CALL f_env_add_defaults(ep_env%f_env_id,f_env=f_env,& 02633 ! new_error=new_error,failure=failure) 02634 ! CALL force_env_get(f_env%force_env,qs_env=sub_qs_env,error=new_error) 02635 ! DO i=1,SIZE(local_sm,2) 02636 ! DO ispin=1,ep_env%nspins 02637 ! IF (ASSOCIATED(local_sm(ispin,i)%matrix)) THEN 02638 02639 ! CALL deallocate_matrix(local_sm(ispin,i)%matrix,error=error) 02640 02641 ! END IF 02642 ! END DO 02643 ! END DO 02644 ! DEALLOCATE(local_sm,stat=stat) 02645 ! CPPostconditionNoFail(stat==0,cp_warning_level,routineP,new_error) 02646 ! CALL f_env_rm_defaults(f_env,new_error,ierr) 02647 ! CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 02648 ! END IF 02649 ! END SUBROUTINE ep_give_back_local_sm 02650 02651 ! ***************************************************************************** 02667 02668 ! CHARACTER(len=*), PARAMETER :: routineN = 'ep_full_to_local_sm', & 02669 ! routineP = moduleN//':'//routineN 02670 02671 ! INTEGER :: disp, i, iblock_col, & 02672 ! iblock_row, ii, ip, irep, j, & 02673 ! lcol, lrow, num_pe, stat, & 02674 ! target_p 02675 ! INTEGER, ALLOCATABLE, DIMENSION(:) :: pos_att, rcv_sizes, rdisp, & 02676 ! sdisp, send_sizes 02677 ! INTEGER, DIMENSION(:), POINTER :: atoms, first_col, first_row, & 02678 ! last_col, last_row 02679 ! LOGICAL :: failure 02680 ! REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf 02681 ! REAL(kind=dp), DIMENSION(:, :), POINTER :: m_block 02682 ! TYPE(real_block_node_type), POINTER :: block_node 02683 ! TYPE(replica_env_type), POINTER :: rep_env 02684 02685 ! failure=.FALSE. 02686 ! NULLIFY(m_block) 02687 ! CPPrecondition(ASSOCIATED(full_sm),cp_failure_level,routineP,error,failure) 02688 ! CPPrecondition(ASSOCIATED(local_sm),cp_failure_level,routineP,error,failure) 02689 ! IF (.NOT.failure) THEN 02690 ! CPPrecondition(SIZE(local_sm,1)<=ispin,cp_failure_level,routineP,error,failure) 02691 ! END IF 02692 ! IF (.NOT. failure) THEN 02693 ! rep_env => ep_env%mol_envs 02694 ! num_pe=rep_env%para_env_inter_rep%num_pe 02695 ! ALLOCATE(send_sizes(0:num_pe-1),stat=stat) 02696 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02697 ! send_sizes=0 02698 ! DO iblock_row=1,full_sm%nblock_row 02699 ! block_node => first_block_node(full_sm,iblock_row) 02700 02701 ! DO WHILE (ASSOCIATED(block_node)) 02702 ! CALL get_block_node(block_node=block_node,& 02703 ! block_col=iblock_col, block=m_block) 02704 ! IF (ep_env%at2sub(iblock_col)==ep_env%at2sub(iblock_row)) THEN 02705 ! target_p=rep_env%inter_rep_rank(& 02706 ! rep_env%replica_owner(ep_env%at2sub(iblock_col))) 02707 ! send_sizes(target_p)=send_sizes(target_p)+SIZE(m_block,1)*SIZE(m_block,2)+2 02708 ! END IF 02709 ! block_node => next_block_node(block_node) 02710 ! END DO 02711 02712 ! END DO 02713 02714 ! ALLOCATE(sdisp(0:num_pe), stat=stat) 02715 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02716 ! disp=0 02717 ! sdisp(0)=0 02718 ! DO ip=1,num_pe 02719 ! disp=disp+send_sizes(ip-1) 02720 ! sdisp(ip)=disp 02721 ! END DO 02722 ! ALLOCATE(send_buf(disp),stat=stat) 02723 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02724 02725 ! ALLOCATE(rcv_sizes(0:num_pe-1),stat=stat) 02726 ! CALL mp_alltoall(send_sizes,rcv_sizes,1,num_pe) 02727 02728 ! ALLOCATE(rdisp(0:num_pe),stat=stat) 02729 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02730 ! disp=0 02731 ! rdisp(0)=0 02732 ! DO ip=1,num_pe 02733 ! disp=disp+rcv_sizes(ip-1) 02734 ! rdisp(ip)=disp 02735 ! END DO 02736 ! ALLOCATE(rcv_buf(disp),stat=stat) 02737 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02738 02739 ! ALLOCATE(pos_att(0:num_pe-1),stat=stat) 02740 ! CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02741 ! pos_att=sdisp(0:num_pe-1) 02742 ! DO iblock_row=1,full_sm%nblock_row 02743 ! block_node => first_block_node(full_sm,iblock_row) 02744 02745 ! DO WHILE (ASSOCIATED(block_node)) 02746 ! CALL get_block_node(block_node=block_node,& 02747 ! block_col=iblock_col,& 02748 ! BLOCK=m_block) 02749 ! IF (ep_env%at2sub(iblock_col)==ep_env%at2sub(iblock_row)) THEN 02750 ! target_p=rep_env%inter_rep_rank(& 02751 ! rep_env%replica_owner(ep_env%at2sub(iblock_col))) 02752 ! send_buf(pos_att(target_p)+1)=REAL(iblock_row,dp) 02753 ! send_buf(pos_att(target_p)+2)=REAL(iblock_col,dp) 02754 ! ii=pos_att(target_p)+2 02755 ! DO j=1,SIZE(m_block,2) 02756 ! DO i=1,SIZE(m_block,1) 02757 ! ii=ii+1 02758 ! send_buf(ii)=m_block(i,j) 02759 ! END DO 02760 ! END DO 02761 ! pos_att(target_p)=ii 02762 ! END IF 02763 ! block_node => next_block_node(block_node) 02764 ! END DO 02765 ! END DO 02766 ! CPPostcondition(ALL(sdisp(1:)==pos_att),cp_failure_level,routineP,error,failure) 02767 02768 ! CALL mp_alltoall ( send_buf, send_sizes, sdisp,& 02769 ! rcv_buf, rcv_sizes, rdisp, rep_env%para_env_inter_rep%group ) 02770 02771 ! CALL get_matrix_info(full_sm,& 02772 ! first_row=first_row,last_row=last_row,& 02773 ! first_col=first_col,last_col=last_col) 02774 02775 ! pos_att=rdisp(0:num_pe-1) 02776 ! DO ip=0,num_pe-1 02777 ! DO 02778 ! IF (pos_att(ip)>=rdisp(ip+1)) EXIT 02779 ! iblock_row=INT(rcv_buf(pos_att(ip)+1)) 02780 ! iblock_col=INT(rcv_buf(pos_att(ip)+2)) 02781 !FM PRINT *,"ip",ip," pos_att(ip)",pos_att(ip),& 02782 !FM " rdisp(ip+1)",rdisp(ip+1)," iblock_col",iblock_col,& 02783 !FM " iblock_row",iblock_row 02784 ! irep=ep_env%at2sub(iblock_row) 02785 ! CPAssert(ep_env%at2sub(iblock_col)==irep,cp_failure_level,routineP,error,failure) 02786 ! ! allow for different basis in full vs sub system, i.e. use proj? 02787 ! ALLOCATE(m_block(last_row(iblock_row)-first_row(iblock_row)+1,& 02788 ! last_col(iblock_col)-first_col(iblock_col)+1),stat=stat) 02789 ! ii=pos_att(ip)+2 02790 ! DO j=1,SIZE(m_block,2) 02791 ! DO i=1,SIZE(m_block,1) 02792 ! ii=ii+1 02793 ! m_block(i,j)=rcv_buf(ii) 02794 ! END DO 02795 ! END DO 02796 ! pos_att(ip)=ii 02797 ! atoms => ep_env%sub_proj(irep)%projection%atoms 02798 02799 ! lrow=iblock_row-atoms(1)+1 02800 ! lcol=iblock_col-atoms(1)+1 02801 ! CPAssert(lrow>0,cp_failure_level,routineP,error,failure) 02802 ! CPAssert(lrow<=SIZE(atoms),cp_failure_level,routineP,error,failure) 02803 ! CPAssert(lcol>0,cp_failure_level,routineP,error,failure) 02804 ! CPAssert(lcol<=SIZE(atoms),cp_failure_level,routineP,error,failure) 02805 ! CPAssert(atoms(lrow)==iblock_row,cp_failure_level,routineP,error,failure) 02806 ! CPAssert(atoms(lcol)==iblock_col,cp_failure_level,routineP,error,failure) 02807 02808 ! CALL add_block_node(matrix=local_sm(ispin,& 02809 ! rep_env_local_index(rep_env,irep,error=error))%matrix,& 02810 ! block_row=lrow,& 02811 ! block_col=lcol,& 02812 ! BLOCK=m_block,error=error) 02813 ! DEALLOCATE(m_block,stat=stat) 02814 ! CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 02815 ! END DO 02816 ! END DO 02817 ! END IF 02818 ! END SUBROUTINE ep_full_to_local_sm 02819 02820 END MODULE ep_methods
1.7.3