CP2K 2.4 (Revision 12889)

qs_scf_post_se.f90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------------!
00002 !   CP2K: A general program to perform molecular dynamics simulations         !
00003 !   Copyright (C) 2000 - 2013  CP2K developers group                          !
00004 !-----------------------------------------------------------------------------!
00005 
00006 ! *****************************************************************************
00013 MODULE qs_scf_post_se
00014   USE atomic_kind_types,               ONLY: atomic_kind_type
00015   USE cell_types,                      ONLY: cell_type
00016   USE cp_array_r_utils,                ONLY: cp_1d_r_p_type
00017   USE cp_control_types,                ONLY: dft_control_type
00018   USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_sparse_matrix
00019   USE cp_dbcsr_types,                  ONLY: cp_dbcsr_p_type,&
00020                                              cp_dbcsr_type
00021   USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
00022                                              cp_fm_struct_release,&
00023                                              cp_fm_struct_type
00024   USE cp_fm_types,                     ONLY: cp_fm_create,&
00025                                              cp_fm_get_info,&
00026                                              cp_fm_init_random,&
00027                                              cp_fm_p_type,&
00028                                              cp_fm_release,&
00029                                              cp_fm_to_fm,&
00030                                              cp_fm_type
00031   USE cp_output_handling,              ONLY: cp_p_file,&
00032                                              cp_print_key_finished_output,&
00033                                              cp_print_key_should_output,&
00034                                              cp_print_key_unit_nr
00035   USE cp_para_types,                   ONLY: cp_para_env_type
00036   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
00037                                              cp_subsys_type
00038   USE distribution_1d_types,           ONLY: distribution_1d_type
00039   USE f77_blas
00040   USE input_constants,                 ONLY: ot_precond_full_all,&
00041                                              state_loc_all
00042   USE input_section_types,             ONLY: section_get_ival,&
00043                                              section_vals_get,&
00044                                              section_vals_get_subs_vals,&
00045                                              section_vals_type,&
00046                                              section_vals_val_get
00047   USE kinds,                           ONLY: dp
00048   USE molecular_states,                ONLY: construct_molecular_states
00049   USE molecule_types_new,              ONLY: molecule_type
00050   USE mulliken,                        ONLY: multipoles
00051   USE particle_list_types,             ONLY: particle_list_type
00052   USE particle_types,                  ONLY: particle_type
00053   USE physcon,                         ONLY: evolt
00054   USE population_analyses,             ONLY: lowdin_population_analysis,&
00055                                              mulliken_population_analysis
00056   USE preconditioner_types,            ONLY: preconditioner_type
00057   USE pw_env_types,                    ONLY: pw_env_get,&
00058                                              pw_env_type
00059   USE pw_pool_types,                   ONLY: pw_pool_create_pw,&
00060                                              pw_pool_give_back_pw,&
00061                                              pw_pool_p_type,&
00062                                              pw_pool_type
00063   USE pw_types,                        ONLY: COMPLEXDATA1D,&
00064                                              REALDATA3D,&
00065                                              REALSPACE,&
00066                                              RECIPROCALSPACE,&
00067                                              pw_p_type
00068   USE qs_environment_types,            ONLY: get_qs_env,&
00069                                              qs_environment_type
00070   USE qs_ks_methods,                   ONLY: qs_ks_did_change,&
00071                                              qs_ks_update_qs_env
00072   USE qs_loc_methods,                  ONLY: qs_loc_driver
00073   USE qs_loc_molecules,                ONLY: wfc_to_molecule
00074   USE qs_loc_types,                    ONLY: localized_wfn_control_create,&
00075                                              qs_loc_env_create,&
00076                                              qs_loc_env_destroy,&
00077                                              qs_loc_env_new_type
00078   USE qs_loc_utils,                    ONLY: qs_loc_env_init,&
00079                                              read_loc_section,&
00080                                              set_loc_centers,&
00081                                              set_loc_wfn_lists
00082   USE qs_mo_methods,                   ONLY: calculate_subspace_eigenvalues
00083   USE qs_mo_types,                     ONLY: get_mo_set,&
00084                                              mo_set_p_type,&
00085                                              set_mo_occupation,&
00086                                              write_mo_set
00087   USE qs_ot_eigensolver,               ONLY: ot_eigensolver
00088   USE qs_resp,                         ONLY: resp_fit
00089   USE qs_rho_types,                    ONLY: qs_rho_type
00090   USE qs_scf_post_gpw,                 ONLY: qs_scf_post_occ_cubes,&
00091                                              qs_scf_post_unocc_cubes
00092   USE qs_scf_types,                    ONLY: qs_scf_env_type
00093   USE scf_control_types,               ONLY: scf_control_type
00094   USE scp_environment_types,           ONLY: get_scp_env,&
00095                                              scp_environment_type
00096   USE timings,                         ONLY: timeset,&
00097                                              timestop
00098 #include "cp_common_uses.h"
00099 
00100   IMPLICIT NONE
00101   PRIVATE
00102 
00103   ! Global parameters
00104   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_se'
00105   PUBLIC :: scf_post_calculation_se
00106 
00107 CONTAINS
00108 
00109 ! *****************************************************************************
00126   SUBROUTINE scf_post_calculation_se(dft_section, scf_env,qs_env,error)
00127 
00128     TYPE(section_vals_type), POINTER         :: dft_section
00129     TYPE(qs_scf_env_type), POINTER           :: scf_env
00130     TYPE(qs_environment_type), POINTER       :: qs_env
00131     TYPE(cp_error_type), INTENT(inout)       :: error
00132 
00133     CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_se', 
00134       routineP = moduleN//':'//routineN
00135 
00136     INTEGER :: handle, homo, ispin, n, n_mo(2), nao, nelectron, nhomo, nlumo, 
00137       nlumos, nmo, nmoloc(2), nmos, output_unit
00138     INTEGER, DIMENSION(:), POINTER           :: marked_states
00139     LOGICAL :: do_mo_cubes, do_wannier_cubes, explicit, failure, ionode, 
00140       my_localized_wfn, p_loc, p_molstates, p_wanstates, uniform_occ
00141     REAL(KIND=dp)                            :: maxocc
00142     REAL(KIND=dp), DIMENSION(:), POINTER     :: mo_eigenvalues
00143     REAL(KIND=dp), DIMENSION(:, :), POINTER  :: scenter
00144     TYPE(cp_1d_r_p_type), DIMENSION(:), 
00145       POINTER                                :: unoccupied_evals
00146     TYPE(cp_dbcsr_p_type), DIMENSION(:), 
00147       POINTER                                :: ks_rmpv, matrix_s
00148     TYPE(cp_dbcsr_type), POINTER             :: mo_coeff_deriv
00149     TYPE(cp_fm_p_type), DIMENSION(:), 
00150       POINTER                                :: unoccupied_orbs
00151     TYPE(cp_fm_struct_type), POINTER         :: fm_struct_tmp
00152     TYPE(cp_fm_type), POINTER                :: mo_coeff, mo_localized
00153     TYPE(cp_logger_type), POINTER            :: logger
00154     TYPE(cp_para_env_type), POINTER          :: para_env
00155     TYPE(cp_subsys_type), POINTER            :: subsys
00156     TYPE(dft_control_type), POINTER          :: dft_control
00157     TYPE(mo_set_p_type), DIMENSION(:), 
00158       POINTER                                :: mos
00159     TYPE(molecule_type), POINTER             :: molecule_set( : )
00160     TYPE(particle_list_type), POINTER        :: particles
00161     TYPE(preconditioner_type), POINTER       :: local_preconditioner
00162     TYPE(pw_env_type), POINTER               :: pw_env
00163     TYPE(pw_p_type)                          :: wf_g, wf_r
00164     TYPE(pw_pool_p_type), DIMENSION(:), 
00165       POINTER                                :: pw_pools
00166     TYPE(pw_pool_type), POINTER              :: auxbas_pw_pool
00167     TYPE(qs_loc_env_new_type), POINTER       :: qs_loc_env
00168     TYPE(qs_rho_type), POINTER               :: rho
00169     TYPE(scf_control_type), POINTER          :: scf_control
00170     TYPE(section_vals_type), POINTER         :: input, loc_print_section, 
00171                                                 localize_section, print_key
00172 
00173     CALL timeset(routineN,handle)
00174 
00175     ! Writes the data that is already available in qs_env
00176     para_env=>qs_env%para_env
00177     CALL write_available_results(qs_env,scf_env,error)
00178 
00179     failure=.FALSE.
00180     my_localized_wfn = .FALSE.
00181     NULLIFY(dft_control, pw_env, auxbas_pw_pool, pw_pools, mos, rho, &
00182          mo_coeff, ks_rmpv, matrix_s, qs_loc_env, scf_control, &
00183          unoccupied_orbs, mo_eigenvalues, unoccupied_evals, fm_struct_tmp, &
00184          molecule_set, subsys, particles, input, print_key)
00185 
00186     logger => cp_error_get_logger(error)
00187     output_unit= cp_logger_get_default_io_unit(logger)
00188 
00189     CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure)
00190     CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure)
00191     CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure)
00192     ! Here we start with data that needs a postprocessing...
00193     IF (.NOT. failure) THEN
00194        CALL get_qs_env(qs_env,dft_control=dft_control,molecule_set=molecule_set,&
00195             mos=mos,rho=rho,matrix_ks=ks_rmpv,scf_control=scf_control,matrix_s=matrix_s,&
00196             input=input, subsys=subsys, pw_env=pw_env, error=error)
00197        CALL cp_subsys_get(subsys,particles=particles,error=error)
00198 
00199        ! Compute Atomic Charges
00200        CALL qs_scf_post_charges(input, logger, qs_env, rho, matrix_s, mos, para_env, error)
00201 
00202        ! Moments of charge distribution
00203        CALL qs_scf_post_moments(input, logger, qs_env, output_unit, error)
00204 
00205        ! Determine if we need to computer properties using the localized centers
00206        localize_section     => section_vals_get_subs_vals(dft_section,"LOCALIZE",error=error)
00207        loc_print_section => section_vals_get_subs_vals(localize_section,"PRINT",error=error)
00208        print_key => section_vals_get_subs_vals(loc_print_section,"MOLECULAR_DIPOLES",error=error)
00209        p_loc=BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00210        print_key => section_vals_get_subs_vals(loc_print_section,"TOTAL_DIPOLE",error=error)
00211        p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00212        print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_CENTERS",error=error)
00213        p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00214        print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_SPREADS",error=error)
00215        p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00216        print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_CUBES",error=error)
00217        p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00218        CALL section_vals_val_get(localize_section,"_SECTION_PARAMETERS_",l_val=explicit,error=error)
00219        p_loc=p_loc.OR.explicit
00220 
00221        ! Print Molecular States
00222        print_key => section_vals_get_subs_vals(loc_print_section,"MOLECULAR_STATES",error=error)
00223        p_molstates=BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00224        p_loc=p_loc.OR.p_molstates
00225        ! Print Wannier States
00226        print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_STATES",error=error)
00227        p_wanstates=BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)
00228        p_loc=p_loc.OR.p_wanstates
00229 
00230        ! Check for cubes MOs and WANNIERS
00231        do_mo_cubes=BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%MO_CUBES",&
00232                    error=error),cp_p_file)
00233        do_wannier_cubes=BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,&
00234                         "WANNIER_CUBES",error=error),cp_p_file)
00235        nlumo=section_get_ival(dft_section,"PRINT%MO_CUBES%NLUMO",error=error)
00236        nhomo=section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO",error=error)
00237 
00238        ! Setup the grids needed to compute a wavefunction given a vector..
00239        IF (((do_mo_cubes.OR.do_wannier_cubes).AND.(nlumo/=0.OR.nhomo/=0.)).OR.p_loc) THEN
00240           CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools,error=error)
00241           CALL pw_pool_create_pw(auxbas_pw_pool,wf_r%pw, use_data = REALDATA3D,&
00242                in_space = REALSPACE, error=error)
00243           CALL pw_pool_create_pw(auxbas_pw_pool,wf_g%pw, use_data = COMPLEXDATA1D,&
00244                in_space = RECIPROCALSPACE, error=error)
00245        END IF
00246 
00247 
00248        ! Makes the MOs eigenstates, computes eigenvalues, write cubes
00249        IF (do_mo_cubes .AND. (nhomo/=0)) THEN
00250           IF (qs_env%dft_control%restricted) THEN
00251              !
00252              ! The issue is that one state is different from the others,
00253              ! so we certainly can't rotate all N+1 alpha states among themselves
00254              ! furthermore, alpha and beta orbitals do have the same shape,
00255              ! but not the same v_xc, therefore these seems to be no reason that
00256              ! they should have the same eigenvalues.
00257              ! So it looks like I need a convincing definition of an eigenstate in this case :-)
00258              !
00259              IF (output_unit>0) WRITE(output_unit,*) &
00260                 " Unclear how we define MOs in the restricted case ... skipping"
00261           ELSE
00262              DO ispin=1,dft_control%nspins
00263                 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff, &
00264                                 eigenvalues=mo_eigenvalues,homo=homo,nmo=nmo)
00265                 IF (output_unit>0) WRITE(output_unit,*) " "
00266                 IF (output_unit>0) WRITE(output_unit,*) " Eigenvalues of the occupied subspace spin ",ispin
00267                 IF (homo .NE. nmo) THEN
00268                    IF (output_unit>0) WRITE(output_unit,*)" and ",nmo-homo," added MO eigenvalues"
00269                 ENDIF
00270                 IF (output_unit>0) WRITE(output_unit,*) "---------------------------------------------"
00271                 ! Also rotate the OT derivs, since they are needed for force calculations
00272                 IF (ASSOCIATED(qs_env%mo_derivs)) THEN
00273                    mo_coeff_deriv=>qs_env%mo_derivs(ispin)%matrix
00274                 ELSE
00275                    mo_coeff_deriv=>NULL()
00276                 ENDIF
00277                 CALL calculate_subspace_eigenvalues(mo_coeff,ks_rmpv(ispin)%matrix,mo_eigenvalues, &
00278                   scr=output_unit, ionode=output_unit>0, do_rotation=.TRUE.,&
00279                   co_rotate_dbcsr=mo_coeff_deriv,error=error)
00280 
00281                 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set,smear=scf_control%smear,error=error)
00282                 IF (output_unit>0) WRITE(output_unit,'(T2,A,F12.6)') "Fermi Energy [eV] :",mos(ispin)%mo_set%mu*evolt
00283                 ! Prints the cube files of OCCUPIED ORBITALS
00284                 CALL qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env,&
00285                      mo_coeff, wf_g, wf_r, particles, homo, ispin, error=error)
00286              ENDDO
00287           ENDIF
00288        ENDIF
00289 
00290        ! Initialize the localization environment, needed e.g. for wannier functions and molecular states
00291        IF (p_loc) THEN
00292           CALL qs_loc_env_create(qs_loc_env,error=error)
00293           CALL localized_wfn_control_create(qs_loc_env%localized_wfn_control,error=error)
00294           CALL read_loc_section(qs_loc_env%localized_wfn_control,localize_section,&
00295                qs_loc_env%do_localize, error=error)
00296           IF (qs_loc_env%do_localize) THEN
00297              ! Some setup for MOs to be localized
00298              nmoloc = 0
00299              DO ispin = 1,dft_control%nspins
00300                 CALL get_mo_set(mos(ispin)%mo_set,nmo=n_mo(ispin),nelectron=nelectron,maxocc=maxocc)
00301                 IF(qs_loc_env%localized_wfn_control%set_of_states == state_loc_all) THEN
00302                    nmoloc(ispin) = NINT(nelectron/maxocc)
00303                 ELSE
00304                    nmoloc(ispin) = MIN(qs_loc_env%localized_wfn_control%nloc_states(1) ,n_mo(ispin))
00305                 ENDIF
00306              ENDDO  ! ispin
00307              CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control,nmoloc,n_mo,dft_control%nspins,error=error)
00308              CALL set_loc_centers(qs_loc_env%localized_wfn_control,nmoloc,dft_control%nspins,error=error)
00309              CALL qs_loc_env_init(qs_loc_env,qs_loc_env%localized_wfn_control,qs_env,error=error)
00310           ELSE
00311              ! Let's inform in case the section is not present in the input
00312              CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,&
00313                   "User requested the calculation of the localized wavefunction but the section "//&
00314                   "LOCALIZE (located in: DFT%LOCALIZE) was not specified. Localization "//&
00315                   "will not be performed!"//&
00316 CPSourceFileRef,&
00317                   only_ionode=.TRUE.)
00318           END IF
00319        END IF
00320 
00321        ! Gets localization info for the occupied orbs
00322        !  - Possibly gets wannier functions
00323        !  - Possibly gets molecular states
00324        IF (p_loc) THEN
00325           IF (qs_env%dft_control%restricted) THEN
00326              IF (output_unit>0) WRITE(output_unit,*) &
00327               " Unclear how we define MOs / localization in the restricted case ... skipping"
00328           ELSE
00329              nmos=0
00330              DO ispin=1,dft_control%nspins
00331                 CALL get_mo_set(mo_set=mos(ispin)%mo_set,nmo=nmo,uniform_occupation=uniform_occ)
00332                 CPPrecondition(uniform_occ,cp_failure_level,routineP,error,failure)
00333                 nmos=MAX(nmos,nmo)
00334              END DO
00335              qs_loc_env%tag_mo="HOMO"
00336              DO ispin=1,dft_control%nspins
00337                 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff, &
00338                                 eigenvalues=mo_eigenvalues,homo=homo,nmo=nmo)
00339                 ! Get eigenstates
00340                 IF (section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO",error=error)==0) THEN
00341                    IF (output_unit>0) WRITE(output_unit,*) " "
00342                    IF (output_unit>0) WRITE(output_unit,*) " Eigenvalues of the occupied subspace spin ",ispin
00343                    IF (homo .NE. nmo) THEN
00344                       IF (output_unit>0) WRITE(output_unit,*) " and ",nmo-homo," added MO eigenvalues"
00345                    ENDIF
00346                    IF (output_unit>0) WRITE(output_unit,*) "---------------------------------------------"
00347                    ! Also rotate the OT derivs, since they are needed for force calculations
00348                    IF (ASSOCIATED(qs_env%mo_derivs)) THEN
00349                       mo_coeff_deriv=>qs_env%mo_derivs(ispin)%matrix
00350                    ELSE
00351                       mo_coeff_deriv=>NULL()
00352                    ENDIF
00353                    CALL calculate_subspace_eigenvalues(mo_coeff,ks_rmpv(ispin)%matrix, &
00354                         evals_arg=mo_eigenvalues, scr=output_unit, &
00355                         ionode=output_unit>0, do_rotation=.TRUE.,error=error)
00356                 END IF
00357                 NULLIFY(mo_localized)
00358                 CALL cp_fm_create(mo_localized,mo_coeff%matrix_struct,error=error)
00359                 CALL cp_fm_to_fm(mo_coeff,mo_localized,error=error)
00360                 IF (qs_loc_env%do_localize) THEN
00361                    ! Do the Real localization..
00362                    IF (output_unit>0) WRITE(output_unit,"(/,T2,A,I3)")&
00363                         "LOCALIZATION| Computing localization properties "//&
00364                         "for OCCUPIED ORBITALS. Spin:",ispin
00365 
00366                    scenter => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
00367                    CALL qs_loc_driver(qs_env,qs_loc_env,localize_section,loc_print_section,&
00368                         myspin=ispin,ext_mo_coeff=mo_localized,error=error)
00369 
00370                    ! maps wfc to molecules, and compute the molecular dipoles if required
00371                    IF (BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,&
00372                        "MOLECULAR_DIPOLES",error=error),cp_p_file) .OR. p_molstates) THEN
00373                       CALL wfc_to_molecule(qs_env, qs_loc_env, loc_print_section, scenter, molecule_set,&
00374                            dft_control%nspins, error)
00375                    END IF
00376 
00377                    ! Compute the molecular states
00378                    IF ( p_molstates ) THEN
00379                       CALL construct_molecular_states(molecule_set, mo_localized, mo_coeff, mo_eigenvalues,&
00380                            ks_rmpv(ispin)%matrix, matrix_s(1)%matrix, qs_env,  wf_r, wf_g,&
00381                            loc_print_section=loc_print_section, particles=particles, tag="HOMO",&
00382                            marked_states=marked_states, error=error)
00383                       IF(ASSOCIATED(marked_states)) DEALLOCATE(marked_states)
00384                    ENDIF
00385                    IF ( p_wanstates ) THEN
00386                       IF(ASSOCIATED(marked_states)) DEALLOCATE(marked_states)
00387                    ENDIF
00388                 END IF
00389                 CALL cp_fm_release(mo_localized,error=error)
00390              ENDDO
00391              ! Print Total Dipole if the localization has been performed
00392              IF (qs_loc_env%do_localize) THEN
00393                 CALL qs_scf_post_loc_dip(input, dft_control, qs_loc_env, logger, qs_env, error)
00394              END IF
00395           ENDIF
00396        ENDIF
00397 
00398        ! Gets the lumos, and eigenvalues for the lumos
00399        IF (do_mo_cubes .AND. nlumo .NE. 0 ) THEN
00400           ALLOCATE(unoccupied_orbs(dft_control%nspins))
00401           ALLOCATE(unoccupied_evals(dft_control%nspins))
00402           DO ispin=1,dft_control%nspins
00403              NULLIFY(unoccupied_orbs(ispin)%matrix)
00404              NULLIFY(unoccupied_evals(ispin)%array)
00405              ! Always write eigenvalues
00406              IF (output_unit>0) WRITE(output_unit,*) " "
00407              IF (output_unit>0) WRITE(output_unit,*) " Lowest Eigenvalues of the unoccupied subspace spin ",ispin
00408              IF (output_unit>0) WRITE(output_unit,FMT='(1X,A)')&
00409                                       "-----------------------------------------------------"
00410              CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff,homo=homo,nao=nao,nmo=nmo)
00411              CALL cp_fm_get_info(mo_coeff, nrow_global=n,error=error)
00412              nlumos=MAX(1,MIN(nlumo,nao-homo))
00413              IF (nlumo==-1) nlumos=nao-homo
00414              ALLOCATE(unoccupied_evals(ispin)%array(nlumos))
00415              CALL cp_fm_struct_create(fm_struct_tmp,para_env=qs_env%para_env,context=qs_env%blacs_env, &
00416                                                     nrow_global=n,ncol_global=nlumos,error=error)
00417              CALL cp_fm_create(unoccupied_orbs(ispin)%matrix, fm_struct_tmp,name="lumos",error=error)
00418              CALL cp_fm_struct_release(fm_struct_tmp,error=error)
00419              CALL cp_fm_init_random(unoccupied_orbs(ispin)%matrix,nlumos,error=error)
00420 
00421              ! the full_all preconditioner makes not much sense for lumos search
00422              NULLIFY(local_preconditioner)
00423              IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
00424                 local_preconditioner=>scf_env%ot_preconditioner(1)%preconditioner
00425                 ! this one can for sure not be right (as it has to match a given C0)
00426                 IF (local_preconditioner%in_use == ot_precond_full_all) THEN
00427                     NULLIFY(local_preconditioner)
00428                 ENDIF
00429              ENDIF
00430 
00431              CALL ot_eigensolver(matrix_h=ks_rmpv(ispin)%matrix,&
00432                   matrix_s=matrix_s(1)%matrix, &
00433                   matrix_c_fm=unoccupied_orbs(ispin)%matrix, &
00434                   matrix_orthogonal_space_fm=mo_coeff, &
00435                   eps_gradient=scf_control%eps_lumos, &
00436                   preconditioner=local_preconditioner, &
00437                   iter_max=scf_control%max_iter_lumos,&
00438                   size_ortho_space=homo,error=error)
00439 
00440              CALL calculate_subspace_eigenvalues(unoccupied_orbs(ispin)%matrix,ks_rmpv(ispin)%matrix,&
00441                   unoccupied_evals(ispin)%array, scr=output_unit, &
00442                   ionode=output_unit>0,error=error)
00443 
00444              ! Prints the cube files of UNOCCUPIED ORBITALS
00445              CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env,&
00446                   unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, error)
00447           ENDDO
00448 
00449        ENDIF
00450 
00451 
00452        ! Deallocate grids needed to compute wavefunctions
00453        IF ( ( ( do_mo_cubes .OR. do_wannier_cubes ).AND. (nlumo /= 0 .OR. nhomo /= 0 )) .OR. p_loc ) THEN
00454           CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_r%pw, error=error)
00455           CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_g%pw, error=error)
00456        END IF
00457 
00458        ! Destroy the localization environment
00459        IF (p_loc) THEN
00460           CALL qs_loc_env_destroy(qs_loc_env, error=error)
00461        ENDIF
00462 
00463        ! Compute the optical conductivity (needs homos and lumos)
00464        CALL qs_scf_post_optc(input, dft_section, dft_control, logger, qs_env,&
00465             unoccupied_orbs, unoccupied_evals, output_unit, error)
00466 
00467        ! This is just a deallocation for printing MO_CUBES
00468        IF (do_mo_cubes .AND. nlumo /= 0 ) THEN
00469           DO ispin=1,dft_control%nspins
00470              DEALLOCATE(unoccupied_evals(ispin)%array)
00471              CALL cp_fm_release(unoccupied_orbs(ispin)%matrix,error=error)
00472           ENDDO
00473           DEALLOCATE(unoccupied_evals)
00474           DEALLOCATE(unoccupied_orbs)
00475        ENDIF
00476 
00477        ! Print coherent X-ray diffraction spectrum
00478        CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit, error)
00479 
00480        ! Calculation of Electric Field Gradients
00481        CALL qs_scf_post_efg(input, logger, error)
00482 
00483        ! Calculation of ET
00484        CALL qs_scf_post_et(input, qs_env, dft_control, error)
00485 
00486        ! Calculation of EPR Hyperfine Coupling Tensors
00487        CALL qs_scf_post_epr(input, logger, qs_env, error)
00488 
00489        ! Calculation of properties needed for BASIS_MOLOPT optimizations
00490        CALL qs_scf_post_molopt(input, logger, qs_env, error)
00491 
00492     END IF
00493 
00494     CALL timestop(handle)
00495 
00496   END SUBROUTINE scf_post_calculation_se
00497 
00498 ! *****************************************************************************
00504   SUBROUTINE qs_scf_post_efg(input, logger, error)
00505     TYPE(section_vals_type), POINTER         :: input
00506     TYPE(cp_logger_type), POINTER            :: logger
00507     TYPE(cp_error_type), INTENT(inout)       :: error
00508 
00509     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_efg', 
00510       routineP = moduleN//':'//routineN
00511 
00512     TYPE(section_vals_type), POINTER         :: print_key
00513 
00514     print_key => section_vals_get_subs_vals(section_vals=input,&
00515                    subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT",&
00516                    error=error)
00517     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),&
00518               cp_p_file)) THEN
00519        CALL cp_unimplemented_error(fromWhere=routineP, &
00520             message="ELECTRIC_FIELD_GRADIENT not implemented for Semi-Empirical calculations!!", &
00521             error=error, error_level=cp_warning_level)
00522     END IF
00523 
00524   END SUBROUTINE qs_scf_post_efg
00525 
00526 ! *****************************************************************************
00532   SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit, error)
00533     TYPE(section_vals_type), POINTER         :: input
00534     TYPE(cp_logger_type), POINTER            :: logger
00535     TYPE(qs_environment_type), POINTER       :: qs_env
00536     INTEGER, INTENT(IN)                      :: output_unit
00537     TYPE(cp_error_type), INTENT(inout)       :: error
00538 
00539     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_moments', 
00540       routineP = moduleN//':'//routineN
00541 
00542     TYPE(section_vals_type), POINTER         :: print_key
00543 
00544     print_key => section_vals_get_subs_vals(section_vals=input,&
00545          subsection_name="DFT%PRINT%MOMENTS",error=error)
00546 
00547     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN
00548        CALL cp_unimplemented_error(fromWhere=routineP, &
00549             message="ELECTRIC MOMENTS not implemented for Semi-Empirical calculations!!", &
00550             error=error, error_level=cp_warning_level)
00551     END IF
00552 
00553   END SUBROUTINE qs_scf_post_moments
00554 
00555 ! *****************************************************************************
00561   SUBROUTINE qs_scf_post_charges(input, logger, qs_env, rho, matrix_s, mos, para_env, error)
00562     TYPE(section_vals_type), POINTER         :: input
00563     TYPE(cp_logger_type), POINTER            :: logger
00564     TYPE(qs_environment_type), POINTER       :: qs_env
00565     TYPE(qs_rho_type), POINTER               :: rho
00566     TYPE(cp_dbcsr_p_type), DIMENSION(:), 
00567       POINTER                                :: matrix_s
00568     TYPE(mo_set_p_type), DIMENSION(:), 
00569       POINTER                                :: mos
00570     TYPE(cp_para_env_type), POINTER          :: para_env
00571     TYPE(cp_error_type), INTENT(inout)       :: error
00572 
00573     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_charges', 
00574       routineP = moduleN//':'//routineN
00575 
00576     INTEGER                                  :: print_level, unit_nr
00577     LOGICAL                                  :: print_it, scp_nddo
00578     TYPE(atomic_kind_type), DIMENSION(:), 
00579       POINTER                                :: atomic_kind_set
00580     TYPE(cp_dbcsr_type), POINTER             :: pscp
00581     TYPE(dft_control_type), POINTER          :: dft_control
00582     TYPE(distribution_1d_type), POINTER      :: local_particles
00583     TYPE(particle_type), DIMENSION(:), 
00584       POINTER                                :: particle_set
00585     TYPE(scp_environment_type), POINTER      :: scp_env
00586     TYPE(section_vals_type), POINTER         :: print_key
00587 
00588     NULLIFY(particle_set)
00589     CALL get_qs_env(qs_env=qs_env,&
00590                     atomic_kind_set=atomic_kind_set,&
00591                     local_particles=local_particles,&
00592                     particle_set=particle_set,&
00593                     scp_env=scp_env,&
00594                     dft_control=dft_control,&
00595                     error=error)
00596     scp_nddo = dft_control%qs_control%se_control%scp
00597     IF (scp_nddo ) &
00598     CALL get_scp_env(scp_env=scp_env, pscp=pscp, error=error)
00599     ! Compute the mulliken charges
00600     print_key => section_vals_get_subs_vals(input,"DFT%PRINT%MULLIKEN", error=error)
00601     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN
00602        CALL cp_unimplemented_error(fromWhere=routineP, &
00603             message="MULLIKEN charges not validated for Semi-Empirical calculations!!", &
00604             error=error, error_level=cp_warning_level)
00605        unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%MULLIKEN",extension=".mulliken",&
00606                                     middle_name="",log_filename=.FALSE.,error=error)
00607        print_level = 1
00608        CALL section_vals_val_get(print_key,"PRINT_GOP",l_val=print_it,error=error)
00609        IF (print_it) print_level = 2
00610        CALL section_vals_val_get(print_key,"PRINT_ALL",l_val=print_it,error=error)
00611        IF (print_it) print_level = 3
00612        CALL mulliken_population_analysis(qs_env,unit_nr,print_level,error)
00613        IF ( scp_nddo ) &
00614        CALL multipoles(atomic_kind_set,particle_set,&
00615          rho%rho_ao,pscp,matrix_s(1)%matrix,para_env,unit_nr,"DFT%PRINT%MULLIKEN",error=error)
00616 
00617        CALL cp_print_key_finished_output(unit_nr, logger,input,"DFT%PRINT%MULLIKEN",error=error)
00618     END IF
00619 
00620     ! Compute the Lowdin charges
00621     print_key => section_vals_get_subs_vals(input,"DFT%PRINT%LOWDIN", error=error)
00622     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN
00623        CALL cp_unimplemented_error(fromWhere=routineP,&
00624                                    message="Lowdin charges not validated for semi-empirical calculations!", &
00625                                    error=error, error_level=cp_warning_level)
00626        unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%LOWDIN",extension=".lowdin",&
00627                                       log_filename=.FALSE.,error=error)
00628        print_level = 1
00629        CALL section_vals_val_get(print_key,"PRINT_GOP",l_val=print_it,error=error)
00630        IF (print_it) print_level = 2
00631        CALL section_vals_val_get(print_key,"PRINT_ALL",l_val=print_it,error=error)
00632        IF (print_it) print_level = 3
00633        CALL lowdin_population_analysis(qs_env,unit_nr,print_level,error)
00634        CALL cp_print_key_finished_output(unit_nr, logger,input,"DFT%PRINT%LOWDIN", error=error)
00635     END IF
00636 
00637     ! Compute the RESP charges
00638     CALL resp_fit(qs_env,error)
00639 
00640     ! Compute the Density Derived Atomic Point charges with the Bloechl scheme
00641     print_key => section_vals_get_subs_vals(input,"PROPERTIES%FIT_CHARGE", error=error)
00642     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN
00643        CALL cp_unimplemented_error(fromWhere=routineP, &
00644             message="FIT CHARGES not implemented for Semi-Empirical calculations!!", &
00645             error=error, error_level=cp_warning_level)
00646     END IF
00647 
00648   END SUBROUTINE qs_scf_post_charges
00649 
00650 ! *****************************************************************************
00656   SUBROUTINE qs_scf_post_loc_dip(input, dft_control, qs_loc_env, logger, qs_env, error)
00657     TYPE(section_vals_type), POINTER         :: input
00658     TYPE(dft_control_type), POINTER          :: dft_control
00659     TYPE(qs_loc_env_new_type), POINTER       :: qs_loc_env
00660     TYPE(cp_logger_type), POINTER            :: logger
00661     TYPE(qs_environment_type), POINTER       :: qs_env
00662     TYPE(cp_error_type), INTENT(inout)       :: error
00663 
00664     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_loc_dip', 
00665       routineP = moduleN//':'//routineN
00666 
00667     LOGICAL                                  :: failure, first_time
00668 
00669     failure = .FALSE.
00670     IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00671          "DFT%PRINT%LOCALIZATION%TOTAL_DIPOLE",first_time=first_time,error=error),cp_p_file)) THEN
00672        CALL cp_unimplemented_error(fromWhere=routineP, &
00673             message="TOTAL DIPOLE  not implemented for Semi-Empirical calculations!!", &
00674             error=error, error_level=cp_warning_level)
00675     END IF
00676 
00677   END SUBROUTINE qs_scf_post_loc_dip
00678 
00679 ! *****************************************************************************
00685   SUBROUTINE qs_scf_post_optc(input, dft_section, dft_control, logger, qs_env,&
00686        unoccupied_orbs, unoccupied_evals, output_unit, error)
00687     TYPE(section_vals_type), POINTER         :: input, dft_section
00688     TYPE(dft_control_type), POINTER          :: dft_control
00689     TYPE(cp_logger_type), POINTER            :: logger
00690     TYPE(qs_environment_type), POINTER       :: qs_env
00691     TYPE(cp_fm_p_type), DIMENSION(:), 
00692       POINTER                                :: unoccupied_orbs
00693     TYPE(cp_1d_r_p_type), DIMENSION(:), 
00694       POINTER                                :: unoccupied_evals
00695     INTEGER, INTENT(IN)                      :: output_unit
00696     TYPE(cp_error_type), INTENT(inout)       :: error
00697 
00698     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_optc', 
00699       routineP = moduleN//':'//routineN
00700 
00701     IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00702          "DFT%PRINT%OPTICAL_CONDUCTIVITY",error=error),cp_p_file)) THEN
00703        CALL cp_unimplemented_error(fromWhere=routineP, &
00704             message="OPTICAL_CONDUCTIVITY  not implemented for Semi-Empirical calculations!!", &
00705             error=error, error_level=cp_warning_level)
00706     ENDIF
00707   END SUBROUTINE qs_scf_post_optc
00708 
00709 ! *****************************************************************************
00715   SUBROUTINE qs_scf_post_xray(input,dft_section,logger,qs_env,&
00716                               output_unit,error)
00717 
00718     TYPE(section_vals_type), POINTER         :: input, dft_section
00719     TYPE(cp_logger_type), POINTER            :: logger
00720     TYPE(qs_environment_type), POINTER       :: qs_env
00721     INTEGER, INTENT(IN)                      :: output_unit
00722     TYPE(cp_error_type), INTENT(inout)       :: error
00723 
00724     CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_post_xray', 
00725       routineP = moduleN//':'//routineN
00726 
00727     TYPE(section_vals_type), POINTER         :: print_key
00728 
00729     print_key => section_vals_get_subs_vals(section_vals=input,&
00730                                             subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM",&
00731                                             error=error)
00732 
00733     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN
00734        CALL cp_unimplemented_error(fromWhere=routineP, &
00735             message="XRAY_DIFFRACTION_SPECTRUM  not implemented for Semi-Empirical calculations!!", &
00736             error=error, error_level=cp_warning_level)
00737     END IF
00738 
00739   END SUBROUTINE qs_scf_post_xray
00740 
00741 ! *****************************************************************************
00747   SUBROUTINE qs_scf_post_epr(input, logger, qs_env, error)
00748     TYPE(section_vals_type), POINTER         :: input
00749     TYPE(cp_logger_type), POINTER            :: logger
00750     TYPE(qs_environment_type), POINTER       :: qs_env
00751     TYPE(cp_error_type), INTENT(inout)       :: error
00752 
00753     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_epr', 
00754       routineP = moduleN//':'//routineN
00755 
00756     TYPE(section_vals_type), POINTER         :: print_key
00757 
00758     print_key => section_vals_get_subs_vals(section_vals=input,&
00759                    subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR",&
00760                    error=error)
00761     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),&
00762               cp_p_file)) THEN
00763        CALL cp_unimplemented_error(fromWhere=routineP, &
00764             message="HYPERFINE_COUPLING_TENSOR  not implemented for Semi-Empirical calculations!!", &
00765             error=error, error_level=cp_warning_level)
00766     END IF
00767 
00768   END SUBROUTINE qs_scf_post_epr
00769 
00770 ! *****************************************************************************
00780   SUBROUTINE qs_scf_post_molopt(input, logger, qs_env, error)
00781     TYPE(section_vals_type), POINTER         :: input
00782     TYPE(cp_logger_type), POINTER            :: logger
00783     TYPE(qs_environment_type), POINTER       :: qs_env
00784     TYPE(cp_error_type), INTENT(inout)       :: error
00785 
00786     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_molopt', 
00787       routineP = moduleN//':'//routineN
00788 
00789     TYPE(section_vals_type), POINTER         :: print_key
00790 
00791     print_key => section_vals_get_subs_vals(section_vals=input,&
00792                    subsection_name="DFT%PRINT%BASIS_MOLOPT_QUANTITIES",&
00793                    error=error)
00794     IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),&
00795               cp_p_file)) THEN
00796        CALL cp_unimplemented_error(fromWhere=routineP, &
00797             message="BASIS_MOLOPT_QUANTITIES  not implemented for Semi-Empirical calculations!!", &
00798             error=error, error_level=cp_warning_level)
00799     END IF
00800 
00801   END SUBROUTINE qs_scf_post_molopt
00802 
00803 ! *****************************************************************************
00809   SUBROUTINE qs_scf_post_et(input, qs_env, dft_control, error)
00810     TYPE(section_vals_type), POINTER         :: input
00811     TYPE(qs_environment_type), POINTER       :: qs_env
00812     TYPE(dft_control_type), POINTER          :: dft_control
00813     TYPE(cp_error_type), INTENT(inout)       :: error
00814 
00815     CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_et', 
00816       routineP = moduleN//':'//routineN
00817 
00818     LOGICAL                                  :: do_et
00819     TYPE(section_vals_type), POINTER         :: et_section
00820 
00821     do_et=.FALSE.
00822     et_section =>  section_vals_get_subs_vals(input,"PROPERTIES%ET_COUPLING",&
00823                                                   error=error)
00824     CALL section_vals_get(et_section,explicit=do_et,error=error)
00825     IF(do_et)THEN
00826        CALL cp_unimplemented_error(fromWhere=routineP, &
00827             message="ET_COUPLING  not implemented for Semi-Empirical calculations!!", &
00828             error=error, error_level=cp_warning_level)
00829     END IF
00830 
00831   END SUBROUTINE qs_scf_post_et
00832 
00833 ! *****************************************************************************
00839   SUBROUTINE write_available_results(qs_env,scf_env,error)
00840     TYPE(qs_environment_type), POINTER       :: qs_env
00841     TYPE(qs_scf_env_type), POINTER           :: scf_env
00842     TYPE(cp_error_type), INTENT(inout)       :: error
00843 
00844     CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results', 
00845       routineP = moduleN//':'//routineN
00846 
00847     INTEGER                                  :: handle, ispin, iw, output_unit
00848     LOGICAL                                  :: failure
00849     TYPE(atomic_kind_type), DIMENSION(:), 
00850       POINTER                                :: atomic_kind_set
00851     TYPE(cell_type), POINTER                 :: cell
00852     TYPE(cp_dbcsr_p_type), DIMENSION(:), 
00853       POINTER                                :: ks_rmpv
00854     TYPE(cp_logger_type), POINTER            :: logger
00855     TYPE(cp_para_env_type), POINTER          :: para_env
00856     TYPE(cp_subsys_type), POINTER            :: subsys
00857     TYPE(dft_control_type), POINTER          :: dft_control
00858     TYPE(mo_set_p_type), DIMENSION(:), 
00859       POINTER                                :: mos
00860     TYPE(molecule_type), POINTER             :: molecule_set( : )
00861     TYPE(particle_list_type), POINTER        :: particles
00862     TYPE(particle_type), DIMENSION(:), 
00863       POINTER                                :: particle_set
00864     TYPE(qs_rho_type), POINTER               :: rho
00865     TYPE(scf_control_type), POINTER          :: scf_control
00866     TYPE(section_vals_type), POINTER         :: dft_section, input
00867 
00868     CALL timeset(routineN,handle)
00869     failure=.FALSE.
00870     NULLIFY(cell, dft_control, mos, atomic_kind_set, particle_set, rho, &
00871          ks_rmpv, scf_control, dft_section, molecule_set, input, &
00872          particles, subsys)
00873     para_env=>qs_env%para_env
00874     logger => cp_error_get_logger(error)
00875     output_unit= cp_logger_get_default_io_unit(logger)
00876 
00877     CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure)
00878     IF (.NOT. failure) THEN
00879        CALL get_qs_env(qs_env,dft_control=dft_control,molecule_set=molecule_set, &
00880             mos=mos,atomic_kind_set=atomic_kind_set,particle_set=particle_set,&
00881             rho=rho,matrix_ks=ks_rmpv,scf_control=scf_control,&
00882             input=input,cell=cell,subsys=subsys,error=error)
00883        CALL cp_subsys_get(subsys,particles=particles,error=error)
00884 
00885        ! *** if the dft_section tells you to do so, write last wavefunction to screen
00886        dft_section => section_vals_get_subs_vals(input,"DFT",error=error)
00887        IF (dft_control%nspins == 2) THEN
00888          CALL write_mo_set(mos(1)%mo_set,atomic_kind_set,particle_set,4,&
00889                            dft_section,spin="ALPHA",last=.TRUE.,error=error)
00890          CALL write_mo_set(mos(2)%mo_set,atomic_kind_set,particle_set,4,&
00891                            dft_section,spin="BETA",last=.TRUE.,error=error)
00892        ELSE
00893          CALL write_mo_set(mos(1)%mo_set,atomic_kind_set,particle_set,4,&
00894                            dft_section,last=.TRUE.,error=error)
00895        END IF
00896 
00897        ! *** at the end of scf print out the projected dos per kind
00898        IF (BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%PDOS",&
00899            error=error),cp_p_file) ) THEN
00900           CALL cp_unimplemented_error(fromWhere=routineP, &
00901                message="PDOS  not implemented for Semi-Empirical calculations!!", &
00902                error=error, error_level=cp_warning_level)
00903        ENDIF
00904 
00905        ! Print the total density (electronic + core charge)
00906        IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00907             "DFT%PRINT%TOT_DENSITY_CUBE", error=error),cp_p_file)) THEN
00908           CALL cp_unimplemented_error(fromWhere=routineP, &
00909                message="TOT_DENSITY_CUBE  not implemented for Semi-Empirical calculations!!", &
00910                error=error, error_level=cp_warning_level)
00911        END IF
00912 
00913        ! Write cube file with electron density
00914        IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00915                  "DFT%PRINT%E_DENSITY_CUBE",error=error),cp_p_file)) THEN
00916           CALL cp_unimplemented_error(fromWhere=routineP, &
00917                message="E_DENSITY_CUBE  not implemented for Semi-Empirical calculations!!", &
00918                error=error, error_level=cp_warning_level)
00919        END IF ! print key
00920 
00921        ! Print the hartree potential
00922        IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00923             "DFT%PRINT%V_HARTREE_CUBE",error=error),cp_p_file)) THEN
00924           CALL cp_unimplemented_error(fromWhere=routineP, &
00925                message="V_HARTREE_CUBE  not implemented for Semi-Empirical calculations!!", &
00926                error=error, error_level=cp_warning_level)
00927        ENDIF
00928 
00929        !    *** write the density matrix ***
00930        IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00931             "DFT%PRINT%AO_MATRICES/DENSITY",error=error),cp_p_file)) THEN
00932           iw = cp_print_key_unit_nr(logger,input,"DFT%PRINT%AO_MATRICES/DENSITY",&
00933                extension=".Log",error=error)
00934           DO ispin=1,dft_control%nspins
00935              CALL cp_dbcsr_write_sparse_matrix(rho%rho_ao(ispin)%matrix,4,6,qs_env,&
00936                   para_env,output_unit=iw,error=error)
00937           END DO
00938           CALL cp_print_key_finished_output(iw,logger,input,&
00939                "DFT%PRINT%AO_MATRICES/DENSITY", error=error)
00940        END IF
00941 
00942        !    **** the Kohn-Sham matrix itself
00943        IF (BTEST(cp_print_key_should_output(logger%iter_info,input,&
00944             "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",error=error),cp_p_file)) THEN
00945           CALL qs_ks_update_qs_env(qs_env%ks_env,qs_env=qs_env,&
00946                error=error,calculate_forces=.FALSE.,just_energy=.FALSE.)
00947           CALL qs_ks_did_change(qs_env%ks_env,rho_changed=.TRUE.,error=error)
00948           iw = cp_print_key_unit_nr(logger,input,"DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",&
00949                extension=".Log",error=error)
00950           CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(1)%matrix,4,6,Qs_env,para_env,output_unit=iw,error=error)
00951           CALL cp_print_key_finished_output(iw,logger,input,&
00952                "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", error=error)
00953        END IF
00954     END IF
00955     CALL timestop(handle)
00956   END SUBROUTINE write_available_results
00957 
00958 END MODULE qs_scf_post_se