|
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 ! ***************************************************************************** 00010 MODULE molecular_states 00011 USE atomic_kind_types, ONLY: atomic_kind_type 00012 USE bibliography, ONLY: Hunt2003,& 00013 cite_reference 00014 USE cell_types, ONLY: cell_type 00015 USE cp_control_types, ONLY: dft_control_type 00016 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply 00017 USE cp_dbcsr_types, ONLY: cp_dbcsr_type 00018 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm 00019 USE cp_fm_diag, ONLY: choose_eigv_solver 00020 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00021 cp_fm_struct_release,& 00022 cp_fm_struct_type 00023 USE cp_fm_types, ONLY: cp_fm_create,& 00024 cp_fm_get_element,& 00025 cp_fm_get_info,& 00026 cp_fm_release,& 00027 cp_fm_to_fm,& 00028 cp_fm_type 00029 USE cp_output_handling, ONLY: cp_p_file,& 00030 cp_print_key_finished_output,& 00031 cp_print_key_should_output,& 00032 cp_print_key_unit_nr 00033 USE cp_para_types, ONLY: cp_para_env_type 00034 USE f77_blas 00035 USE input_section_types, ONLY: section_get_ivals,& 00036 section_vals_type,& 00037 section_vals_val_get 00038 USE kinds, ONLY: default_path_length,& 00039 dp,& 00040 dp_size,& 00041 int_size 00042 USE message_passing, ONLY: mp_bcast,& 00043 mp_maxloc 00044 USE molecule_types_new, ONLY: molecule_type 00045 USE particle_list_types, ONLY: particle_list_type 00046 USE particle_types, ONLY: particle_type 00047 USE pw_env_types, ONLY: pw_env_type 00048 USE pw_types, ONLY: pw_p_type 00049 USE qs_collocate_density, ONLY: calculate_wavefunction 00050 USE qs_environment_types, ONLY: get_qs_env,& 00051 qs_environment_type 00052 USE realspace_grid_cube, ONLY: pw_to_cube 00053 USE termination, ONLY: stop_memory 00054 USE timings, ONLY: timeset,& 00055 timestop 00056 #include "cp_common_uses.h" 00057 00058 IMPLICIT NONE 00059 00060 PRIVATE 00061 00062 ! *** Global parameters *** 00063 00064 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_states' 00065 00066 LOGICAL, PARAMETER, PRIVATE :: debug_this_module=.FALSE. 00067 00068 ! *** Public subroutines *** 00069 00070 PUBLIC :: construct_molecular_states 00071 00072 CONTAINS 00073 00074 ! ***************************************************************************** 00077 SUBROUTINE construct_molecular_states(molecule_set, mo_localized, & 00078 mo_coeff, mo_eigenvalues, Hks, matrix_S, qs_env, wf_r, wf_g,& 00079 loc_print_section, particles, tag, marked_states,error) 00080 00081 TYPE(molecule_type), DIMENSION(:), 00082 POINTER :: molecule_set 00083 TYPE(cp_fm_type), POINTER :: mo_localized, mo_coeff 00084 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues 00085 TYPE(cp_dbcsr_type), POINTER :: Hks, matrix_S 00086 TYPE(qs_environment_type), POINTER :: qs_env 00087 TYPE(pw_p_type), INTENT(INOUT) :: wf_r, wf_g 00088 TYPE(section_vals_type), POINTER :: loc_print_section 00089 TYPE(particle_list_type), POINTER :: particles 00090 CHARACTER(LEN=*), INTENT(IN) :: tag 00091 INTEGER, DIMENSION(:), POINTER :: marked_states 00092 TYPE(cp_error_type), INTENT(INOUT) :: error 00093 00094 CHARACTER(len=*), PARAMETER :: routineN = 'construct_molecular_states', 00095 routineP = moduleN//':'//routineN 00096 00097 CHARACTER(LEN=80) :: title 00098 CHARACTER(LEN=default_path_length) :: filename 00099 INTEGER :: handle, i, imol, iproc, isos, k, n_rep, ncol_global, nproc, 00100 nrow_global, ns, output_unit, unit_nr, unit_report 00101 INTEGER, DIMENSION(:), POINTER :: ind, mark_list 00102 INTEGER, DIMENSION(:, :), POINTER :: mark_states 00103 INTEGER, POINTER :: nstates(:), states(:) 00104 LOGICAL :: explicit, failure 00105 REAL(KIND=dp) :: tmp 00106 REAL(KIND=dp), ALLOCATABLE :: evals( : ) 00107 REAL(KIND=dp), DIMENSION(:), POINTER :: eval_range 00108 TYPE(atomic_kind_type), DIMENSION(:), 00109 POINTER :: atomic_kind_set 00110 TYPE(cell_type), POINTER :: cell 00111 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00112 TYPE(cp_fm_type), POINTER :: b, c, d, D_igk, e_vectors, 00113 rot_e_vectors, smo, storage 00114 TYPE(cp_logger_type), POINTER :: logger 00115 TYPE(cp_para_env_type), POINTER :: para_env 00116 TYPE(dft_control_type), POINTER :: dft_control 00117 TYPE(particle_type), DIMENSION(:), 00118 POINTER :: particle_set 00119 TYPE(pw_env_type), POINTER :: pw_env 00120 00121 CALL cite_reference(Hunt2003) 00122 00123 CALL timeset(routineN,handle) 00124 00125 failure = .FALSE. 00126 NULLIFY(logger,mark_states,mark_list) 00127 logger => cp_error_get_logger(error) 00128 !----------------------------------------------------------------------------- 00129 ! 1. 00130 !----------------------------------------------------------------------------- 00131 para_env => qs_env % para_env 00132 nproc = para_env%num_pe 00133 output_unit = cp_logger_get_default_io_unit(logger) 00134 CALL section_vals_val_get(loc_print_section,"MOLECULAR_STATES%CUBE_EVAL_RANGE",& 00135 explicit=explicit,error=error) 00136 IF (explicit) THEN 00137 CALL section_vals_val_get(loc_print_section,"MOLECULAR_STATES%CUBE_EVAL_RANGE",& 00138 r_vals=eval_range,error=error) 00139 ELSE 00140 ALLOCATE(eval_range(2)) 00141 eval_range(1)=-HUGE(0.0_dp) 00142 eval_range(2)=+HUGE(0.0_dp) 00143 ENDIF 00144 CALL section_vals_val_get(loc_print_section,"MOLECULAR_STATES%MARK_STATES", & 00145 n_rep_val=n_rep, error=error) 00146 IF(n_rep.GT.0)THEN 00147 ALLOCATE (mark_states(2,n_rep),STAT=isos) 00148 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00149 "mark_states",int_size*2*n_rep) 00150 IF (.NOT.ASSOCIATED(marked_states)) THEN 00151 ALLOCATE (marked_states(n_rep),stat=isos) 00152 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00153 "mark_states",int_size*n_rep) 00154 END IF 00155 DO i=1,n_rep 00156 CALL section_vals_val_get(loc_print_section,"MOLECULAR_STATES%MARK_STATES",& 00157 i_rep_val=i,i_vals=mark_list, error=error) 00158 mark_states(:,i)=mark_list(:) 00159 END DO 00160 ELSE 00161 NULLIFY(marked_states) 00162 END IF 00163 00164 !----------------------------------------------------------------------------- 00165 !----------------------------------------------------------------------------- 00166 ! 2. 00167 !----------------------------------------------------------------------------- 00168 unit_report=cp_print_key_unit_nr(logger,loc_print_section,"MOLECULAR_STATES",& 00169 extension=".data",middle_name="Molecular_DOS",log_filename=.FALSE.,error=error) 00170 IF (unit_report>0) THEN 00171 WRITE(unit_report,*) SIZE(mo_eigenvalues)," number of states " 00172 DO i=1,SIZE(mo_eigenvalues) 00173 WRITE(unit_report,*) mo_eigenvalues(i) 00174 ENDDO 00175 ENDIF 00176 00177 !----------------------------------------------------------------------------- 00178 !----------------------------------------------------------------------------- 00179 ! 3. 00180 !----------------------------------------------------------------------------- 00181 CALL cp_fm_get_info(mo_localized, & 00182 ncol_global=ncol_global, & 00183 nrow_global=nrow_global ,error=error) 00184 NULLIFY(smo) 00185 CALL cp_fm_create(smo,mo_coeff%matrix_struct,error=error) 00186 CALL cp_dbcsr_sm_fm_multiply(matrix_S,mo_coeff,smo,ncol_global,error=error) 00187 00188 !----------------------------------------------------------------------------- 00189 !----------------------------------------------------------------------------- 00190 ! 4. 00191 !----------------------------------------------------------------------------- 00192 ALLOCATE (nstates(2),STAT=isos) 00193 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00194 "nstates",int_size*2) 00195 00196 CALL cp_fm_create(storage,mo_localized%matrix_struct,name='storage',error=error) 00197 00198 DO imol = 1, SIZE(molecule_set) 00199 IF (ASSOCIATED(molecule_set(imol)%lmi)) THEN 00200 nstates(1) = molecule_set(imol)%lmi%nstates 00201 ELSE 00202 nstates(1) = 0 00203 ENDIF 00204 nstates(2) = para_env%mepos 00205 00206 CALL mp_maxloc(nstates,para_env%group) 00207 00208 IF (nstates(1)==0) CYCLE 00209 NULLIFY(states) 00210 ALLOCATE(states(nstates(1)),STAT=isos) 00211 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00212 "states",int_size*nstates(1)) 00213 states(:) = 0 00214 00215 iproc=nstates(2) 00216 IF(iproc == para_env%mepos) THEN 00217 states(:) = molecule_set(imol)%lmi%states(:) 00218 END IF 00219 !!BCAST from here root = iproc 00220 CALL mp_bcast(states,iproc,para_env%group) 00221 00222 ns = nstates(1) 00223 ind => states( : ) 00224 ALLOCATE ( evals ( ns ), STAT = isos ) 00225 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00226 "evals",dp_size*ns) 00227 00228 NULLIFY(b,c,d,fm_struct_tmp,e_vectors,rot_e_vectors) 00229 00230 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_global, & 00231 ncol_global=ns, & 00232 para_env=mo_localized%matrix_struct%para_env, & 00233 context=mo_localized%matrix_struct%context,error=error) 00234 00235 CALL cp_fm_create(b,fm_struct_tmp, name="b",error=error) 00236 CALL cp_fm_create(c,fm_struct_tmp, name="c",error=error) 00237 CALL cp_fm_create(rot_e_vectors,fm_struct_tmp, name="rot_e_vectors",error=error) 00238 00239 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00240 00241 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ns, ncol_global=ns, & 00242 para_env=mo_localized%matrix_struct%para_env, & 00243 context=mo_localized%matrix_struct%context,error=error) 00244 00245 CALL cp_fm_create(d,fm_struct_tmp, name="d",error=error) 00246 CALL cp_fm_create(e_vectors,fm_struct_tmp, name="e_vectors",error=error) 00247 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00248 00249 DO i=1,ns 00250 CALL cp_fm_to_fm ( mo_localized, b, 1, ind ( i ), i) 00251 END DO 00252 00253 CALL cp_dbcsr_sm_fm_multiply(Hks,b,c,ns,error=error) 00254 00255 CALL cp_fm_gemm('T','N',ns,ns,nrow_global,1.0_dp, & 00256 b,c,0.0_dp,d,error=error) 00257 00258 CALL choose_eigv_solver( d, e_vectors, evals ,error=error) 00259 00260 IF (output_unit>0) WRITE(output_unit,*)"" 00261 IF (output_unit>0) WRITE(output_unit,*)"MOLECULE ",imol 00262 IF (output_unit>0) WRITE(output_unit,*)"NUMBER OF STATES ", ns 00263 IF (output_unit>0) WRITE(output_unit,*)"EIGENVALUES" 00264 IF (output_unit>0) WRITE(output_unit,*)"" 00265 IF (output_unit>0) WRITE(output_unit,*)"ENERGY original MO-index" 00266 00267 DO k=1,ns 00268 IF(ASSOCIATED(mark_states))THEN 00269 DO i=1,n_rep 00270 IF(imol==mark_states(1,i).AND.k==mark_states(2,i))marked_states(i)=ind(k) 00271 END DO 00272 END IF 00273 IF (output_unit>0) WRITE(output_unit,*) evals(k), ind(k) 00274 END DO 00275 IF (unit_report>0) THEN 00276 WRITE(unit_report,*) imol,ns, " imol, number of states" 00277 DO k=1,ns 00278 WRITE(unit_report,*) evals(k) 00279 ENDDO 00280 ENDIF 00281 00282 CALL cp_fm_gemm('N','N',nrow_global,ns,ns,1.0_dp, & 00283 b,e_vectors,0.0_dp,rot_e_vectors,error=error) 00284 00285 DO i=1,ns 00286 CALL cp_fm_to_fm ( rot_e_vectors, storage, 1, i, ind ( i )) 00287 END DO 00288 00289 IF (.FALSE.) THEN ! this is too much data for large systems 00290 ! compute Eq. 6 from P. Hunt et al. (CPL 376, p. 68-74) 00291 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ns, & 00292 ncol_global=ncol_global, & 00293 para_env=mo_localized%matrix_struct%para_env, & 00294 context=mo_localized%matrix_struct%context,error=error) 00295 CALL cp_fm_create(D_igk,fm_struct_tmp,error=error) 00296 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00297 CALL cp_fm_gemm('T','N',ns,ncol_global,nrow_global,1.0_dp, & 00298 rot_e_vectors,smo,0.0_dp,D_igk,error=error) 00299 DO i=1,ns 00300 DO k=1,ncol_global 00301 CALL cp_fm_get_element(D_igk,i,k,tmp) 00302 IF (unit_report>0) THEN 00303 WRITE(unit_report,*) tmp**2 00304 ENDIF 00305 ENDDO 00306 ENDDO 00307 CALL cp_fm_release(D_igk,error=error) 00308 ENDIF 00309 00310 IF ( BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00311 "MOLECULAR_STATES%CUBES",error=error),cp_p_file) ) THEN 00312 00313 CALL get_qs_env(qs_env=qs_env,& 00314 atomic_kind_set=atomic_kind_set,& 00315 cell=cell,& 00316 dft_control=dft_control,& 00317 particle_set=particle_set,& 00318 pw_env=pw_env,error=error) 00319 00320 DO i=1,ns 00321 IF (evals(i)<eval_range(1) .OR. evals(i)>eval_range(2)) CYCLE 00322 00323 CALL calculate_wavefunction(rot_e_vectors,i,wf_r, & 00324 wf_g, atomic_kind_set,cell,dft_control,particle_set, & 00325 pw_env, error = error ) 00326 00327 WRITE(filename,'(a9,I4.4,a1,I5.5,a4)')"MOLECULE_",imol,"_",i,tag 00328 WRITE(title,'(A,I0,A,I0,A,F14.6,A,I0)') "Mol. Eigenstate ",i," of ",ns," E [a.u.] = ",& 00329 evals(i)," Orig. index ",ind(i) 00330 unit_nr=cp_print_key_unit_nr(logger,loc_print_section,"MOLECULAR_STATES%CUBES",& 00331 extension=".cube",middle_name=TRIM(filename),log_filename=.FALSE.,error=error) 00332 CALL pw_to_cube(wf_r%pw,unit_nr,particles=particles,title=title,& 00333 stride=section_get_ivals(loc_print_section,& 00334 "MOLECULAR_STATES%CUBES%STRIDE",error=error), error=error) 00335 CALL cp_print_key_finished_output(unit_nr,logger,loc_print_section,& 00336 "MOLECULAR_STATES%CUBES",error=error) 00337 END DO 00338 ENDIF 00339 00340 DEALLOCATE (evals,STAT=isos) 00341 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"eval") 00342 CALL cp_fm_release ( b ,error=error) 00343 CALL cp_fm_release ( c ,error=error) 00344 CALL cp_fm_release ( d ,error=error) 00345 CALL cp_fm_release ( e_vectors ,error=error) 00346 CALL cp_fm_release ( rot_e_vectors ,error=error) 00347 00348 DEALLOCATE (states,STAT=isos) 00349 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"states") 00350 00351 END DO 00352 CALL cp_fm_release(smo,error=error) 00353 CALL cp_fm_to_fm(storage,mo_localized,error) 00354 CALL cp_fm_release(storage,error) 00355 IF (ASSOCIATED(mark_states))THEN 00356 DEALLOCATE (mark_states,stat=isos) 00357 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"mark_states") 00358 END IF 00359 DEALLOCATE (nstates,STAT=isos) 00360 IF (isos /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"nstates") 00361 CALL cp_print_key_finished_output(unit_report,logger,loc_print_section,& 00362 "MOLECULAR_STATES",error=error) 00363 !------------------------------------------------------------------------------ 00364 00365 IF (.NOT.explicit) THEN 00366 DEALLOCATE(eval_range) 00367 END IF 00368 00369 CALL timestop(handle) 00370 00371 END SUBROUTINE construct_molecular_states 00372 00373 END MODULE molecular_states
1.7.3