CP2K 2.4 (Revision 12889)

qs_linres_epr_utils.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 ! *****************************************************************************
00020 MODULE qs_linres_epr_utils
00021   USE atomic_kind_types,               ONLY: atomic_kind_type
00022   USE cell_types,                      ONLY: cell_type
00023   USE cp_control_types,                ONLY: dft_control_type
00024   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
00025                                              cp_print_key_unit_nr
00026   USE f77_blas
00027   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
00028                                              section_vals_type
00029   USE kinds,                           ONLY: dp
00030   USE mathconstants,                   ONLY: fourpi,&
00031                                              twopi
00032   USE particle_types,                  ONLY: particle_type
00033   USE physcon,                         ONLY: a_fine,&
00034                                              e_gfactor
00035   USE pw_env_types,                    ONLY: pw_env_get,&
00036                                              pw_env_type
00037   USE pw_pool_types,                   ONLY: pw_pool_create_pw,&
00038                                              pw_pool_type
00039   USE pw_types,                        ONLY: COMPLEXDATA1D,&
00040                                              REALDATA3D,&
00041                                              REALSPACE,&
00042                                              RECIPROCALSPACE,&
00043                                              pw_release
00044   USE qs_environment_types,            ONLY: get_qs_env,&
00045                                              qs_environment_type
00046   USE qs_linres_types,                 ONLY: deallocate_nablavks_atom_set,&
00047                                              epr_env_create,&
00048                                              epr_env_type,&
00049                                              init_nablavks_atom_set,&
00050                                              linres_control_type,&
00051                                              nablavks_atom_type,&
00052                                              set_epr_env
00053   USE qs_matrix_pools,                 ONLY: qs_matrix_pools_type
00054   USE qs_mo_types,                     ONLY: mo_set_p_type
00055   USE qs_rho_atom_types,               ONLY: deallocate_rho_atom_set
00056   USE scf_control_types,               ONLY: scf_control_type
00057   USE timings,                         ONLY: timeset,&
00058                                              timestop
00059 #include "cp_common_uses.h"
00060 
00061   IMPLICIT NONE
00062 
00063   PRIVATE
00064 
00065   PUBLIC :: epr_env_cleanup, epr_env_init
00066 
00067   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_utils'
00068 
00069 CONTAINS
00070 
00071 ! *****************************************************************************
00077   SUBROUTINE epr_env_init(epr_env,qs_env,error)
00078     !
00079     TYPE(epr_env_type)                       :: epr_env
00080     TYPE(qs_environment_type), POINTER       :: qs_env
00081     TYPE(cp_error_type), INTENT(INOUT)       :: error
00082 
00083     CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_env_init', 
00084       routineP = moduleN//':'//routineN
00085 
00086     INTEGER                                  :: handle, i_B, idir, ispin, 
00087                                                 istat, n_mo(2), nao, natom, 
00088                                                 nmoloc, nspins, output_unit
00089     LOGICAL                                  :: failure, gapw
00090     TYPE(atomic_kind_type), DIMENSION(:), 
00091       POINTER                                :: atomic_kind_set
00092     TYPE(cell_type), POINTER                 :: cell
00093     TYPE(cp_logger_type), POINTER            :: logger
00094     TYPE(dft_control_type), POINTER          :: dft_control
00095     TYPE(linres_control_type), POINTER       :: linres_control
00096     TYPE(mo_set_p_type), DIMENSION(:), 
00097       POINTER                                :: mos
00098     TYPE(nablavks_atom_type), DIMENSION(:), 
00099       POINTER                                :: nablavks_atom_set
00100     TYPE(particle_type), DIMENSION(:), 
00101       POINTER                                :: particle_set
00102     TYPE(pw_env_type), POINTER               :: pw_env
00103     TYPE(pw_pool_type), POINTER              :: auxbas_pw_pool
00104     TYPE(qs_matrix_pools_type), POINTER      :: mpools
00105     TYPE(scf_control_type), POINTER          :: scf_control
00106     TYPE(section_vals_type), POINTER         :: lr_section
00107 
00108     CALL timeset(routineN,handle)
00109 
00110     failure = .FALSE.
00111 
00112     NULLIFY(atomic_kind_set, cell, dft_control, linres_control, scf_control)
00113     NULLIFY(logger, mos, mpools, particle_set)
00114     NULLIFY(auxbas_pw_pool, pw_env)
00115     NULLIFY(nablavks_atom_set)
00116 
00117     n_mo(1:2) = 0
00118     nao = 0
00119     nmoloc = 0
00120 
00121     logger => cp_error_get_logger(error)
00122     !ionode = logger%para_env%mepos==logger%para_env%source
00123     lr_section => section_vals_get_subs_vals(qs_env%input,"PROPERTIES%LINRES",error=error)
00124 
00125     output_unit = cp_print_key_unit_nr(logger,lr_section,"PRINT%PROGRAM_RUN_INFO",&
00126          &                             extension=".linresLog",error=error)
00127 
00128     IF(epr_env%ref_count /= 0) THEN
00129       CALL epr_env_cleanup(epr_env,error=error)
00130     END IF
00131 
00132     IF(output_unit>0) THEN
00133       WRITE(output_unit,"(/,T20,A,/)") "*** Start EPR g tensor calculation ***"
00134       WRITE(output_unit,"(T10,A,/)") "Initialization of the EPR environment"
00135     ENDIF
00136 
00137     CALL epr_env_create(epr_env,error=error)
00138 
00139 
00140     CALL get_qs_env(qs_env=qs_env,&
00141                     atomic_kind_set=atomic_kind_set,&
00142                     cell=cell,&
00143                     dft_control=dft_control,&
00144                     linres_control=linres_control,&
00145                     mos=mos,&
00146                     mpools=mpools,&
00147                     particle_set=particle_set,&
00148                     pw_env=pw_env,&
00149                     scf_control=scf_control,error=error)
00150     !
00151     ! Check if restat also psi0 should be restarted
00152     !IF(epr_env%restart_epr .AND. scf_control%density_guess/=restart_guess)THEN
00153     !   CALL stop_program(routineN,moduleN,__LINE__,"restart_epr requires density_guess=restart")
00154     !ENDIF
00155     !
00156     ! check that the psi0 are localized and you have all the centers
00157     CPPrecondition(linres_control%localized_psi0,cp_warning_level,routineP,error,failure)
00158     IF(failure .AND. (output_unit>0)) THEN
00159        WRITE(output_unit,'(A)') &
00160             ' To get EPR parameters within PBC you need localized zero order orbitals '
00161     ENDIF
00162     gapw = dft_control%qs_control%gapw
00163     nspins = dft_control%nspins
00164     natom = SIZE(particle_set,1)
00165     !
00166     ! Conversion factors
00167     ! Magical constant twopi/cell%deth just like in NMR shift (basically undo scale_fac in qs_linres_nmr_current.F)
00168     epr_env%g_free_factor = -1.0_dp * e_gfactor
00169     epr_env%g_zke_factor = e_gfactor * ( a_fine )**2
00170     epr_env%g_so_factor = ( a_fine )**2 * ( -1.0_dp * e_gfactor - 1.0_dp ) / 2.0_dp * twopi / cell%deth
00171     epr_env%g_so_factor_gapw =  ( a_fine )**2 * ( -1.0_dp * e_gfactor - 1.0_dp ) / 2.0_dp
00172     ! * 2 because B_ind = 2 * B_beta
00173     epr_env%g_soo_factor = 2.0_dp * fourpi * ( a_fine )**2 * twopi / cell%deth
00174     ! 2 * 2 * 1/4 * e^2 / m * a_0^2 * 2/3 * mu_0 / (omega * 1e-30 )
00175     epr_env%g_soo_chicorr_factor = 2.0 / 3.0_dp * fourpi * ( a_fine )**2 / cell%deth
00176     !
00177     ! If the current density on the grid needs to be stored
00178     CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,error=error)
00179     !
00180     ! Initialize local current density if GAPW calculation
00181     IF(gapw) THEN
00182        CALL init_nablavks_atom_set(nablavks_atom_set,atomic_kind_set,&
00183                    nspins,error=error)
00184        CALL set_epr_env(epr_env=epr_env,&
00185                         nablavks_atom_set=nablavks_atom_set,&
00186             &           error=error)
00187     ENDIF
00188     !
00189     ! Bind
00190     ALLOCATE(epr_env%bind_set(3,3),STAT=istat)
00191     CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00192     DO i_B = 1,3
00193       DO idir = 1,3
00194         NULLIFY(epr_env%bind_set(idir,i_B)%rho)
00195         ALLOCATE(epr_env%bind_set(idir,i_B)%rho,STAT=istat)
00196         CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00197         epr_env%bind_set(idir,i_B)%rho%ref_count = 1
00198         NULLIFY( epr_env%bind_set(idir,i_B)%rho%rho_r)
00199         NULLIFY( epr_env%bind_set(idir,i_B)%rho%rho_g)
00200 
00201         ALLOCATE(epr_env%bind_set(idir,i_B)%rho%rho_r(1),stat=istat)
00202         CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00203         ALLOCATE(epr_env%bind_set(idir,i_B)%rho%rho_g(1),stat=istat)
00204         CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00205         CALL pw_pool_create_pw(auxbas_pw_pool,&
00206              epr_env%bind_set(idir,i_B)%rho%rho_r(1)%pw,&
00207              use_data=REALDATA3D,in_space=REALSPACE,error=error)
00208         CALL pw_pool_create_pw(auxbas_pw_pool,&
00209              epr_env%bind_set(idir,i_B)%rho%rho_g(1)%pw,&
00210              use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,error=error)
00211       END DO
00212     END DO
00213 
00214     ! Nabla_V_ks
00215     ALLOCATE(epr_env%nablavks_set(3,dft_control%nspins),STAT=istat)
00216     CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00217     DO idir = 1,3
00218       DO ispin = 1,nspins
00219         NULLIFY(epr_env%nablavks_set(idir,ispin)%rho)
00220         ALLOCATE(epr_env%nablavks_set(idir,ispin)%rho,STAT=istat)
00221         CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00222         epr_env%nablavks_set(idir,ispin)%rho%ref_count = 1
00223         NULLIFY( epr_env%nablavks_set(idir,ispin)%rho%rho_r)
00224         NULLIFY( epr_env%nablavks_set(idir,ispin)%rho%rho_g)
00225         ALLOCATE(epr_env%nablavks_set(idir,ispin)%rho%rho_r(1),stat=istat)
00226         CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00227         ALLOCATE(epr_env%nablavks_set(idir,ispin)%rho%rho_g(1),stat=istat)
00228         CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00229         CALL pw_pool_create_pw(auxbas_pw_pool,&
00230              epr_env%nablavks_set(idir,ispin)%rho%rho_r(1)%pw,&
00231              use_data=REALDATA3D,in_space=REALSPACE,error=error)
00232         CALL pw_pool_create_pw(auxbas_pw_pool,&
00233              epr_env%nablavks_set(idir,ispin)%rho%rho_g(1)%pw,&
00234              use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,error=error)
00235       END DO
00236     END DO
00237 
00238     ! Initialize the g tensor components
00239     ALLOCATE(epr_env%g_total(3,3),STAT=istat)
00240     CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00241     ALLOCATE(epr_env%g_so(3,3),STAT=istat)
00242     CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00243     ALLOCATE(epr_env%g_soo(3,3),STAT=istat)
00244     CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00245     epr_env%g_total = 0.0_dp
00246     epr_env%g_zke = 0.0_dp
00247     epr_env%g_so = 0.0_dp
00248     epr_env%g_soo = 0.0_dp
00249 
00250     CALL cp_print_key_finished_output(output_unit,logger,lr_section,&
00251          &                            "PRINT%PROGRAM_RUN_INFO",error=error)
00252 
00253     CALL timestop(handle)
00254 
00255   END SUBROUTINE epr_env_init
00256 
00257 ! *****************************************************************************
00263   SUBROUTINE epr_env_cleanup(epr_env,error)
00264 
00265     TYPE(epr_env_type)                       :: epr_env
00266     TYPE(cp_error_type), INTENT(INOUT)       :: error
00267 
00268     CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_env_cleanup', 
00269       routineP = moduleN//':'//routineN
00270 
00271     INTEGER                                  :: i, i_B, idir, ispin, istat
00272     LOGICAL                                  :: failure
00273 
00274     failure=.FALSE.
00275     IF(.NOT. failure) THEN
00276       epr_env%ref_count = epr_env%ref_count - 1
00277       IF(epr_env%ref_count == 0 ) THEN
00278         ! nablavks_set
00279         IF(ASSOCIATED(epr_env%nablavks_set)) THEN
00280            DO ispin = 1,SIZE(epr_env%nablavks_set,2)
00281            DO idir = 1,SIZE(epr_env%nablavks_set,1)
00282               DO i=1,SIZE(epr_env%nablavks_set(idir,ispin)%rho%rho_r)
00283                  CALL pw_release(epr_env%nablavks_set(idir,ispin)%rho%rho_r(i)%pw,error=error)
00284               END DO
00285               DO i=1,SIZE(epr_env%nablavks_set(idir,ispin)%rho%rho_g)
00286                  CALL pw_release(epr_env%nablavks_set(idir,ispin)%rho%rho_g(i)%pw,error=error)
00287               END DO
00288               DEALLOCATE(epr_env%nablavks_set(idir,ispin)%rho%rho_g,stat=istat)
00289               CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00290               DEALLOCATE(epr_env%nablavks_set(idir,ispin)%rho%rho_r,stat=istat)
00291               CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00292               DEALLOCATE(epr_env%nablavks_set(idir,ispin)%rho,STAT=istat)
00293               CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00294            ENDDO
00295            ENDDO
00296            DEALLOCATE(epr_env%nablavks_set,STAT=istat)
00297            CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00298         END IF
00299         ! nablavks_atom_set
00300         IF(ASSOCIATED(epr_env%nablavks_atom_set)) THEN
00301           CALL deallocate_nablavks_atom_set(epr_env%nablavks_atom_set,error=error)
00302         END IF
00303         ! vks_atom_set
00304         IF(ASSOCIATED(epr_env%vks_atom_set)) THEN
00305           CALL deallocate_rho_atom_set(epr_env%vks_atom_set)
00306         END IF
00307         ! bind_set
00308         IF(ASSOCIATED(epr_env%bind_set)) THEN
00309            DO i_B = 1,SIZE(epr_env%bind_set,2)
00310            DO idir = 1,SIZE(epr_env%bind_set,1)
00311               DO i=1,SIZE(epr_env%bind_set(idir,i_B)%rho%rho_r)
00312                  CALL pw_release(epr_env%bind_set(idir,i_B)%rho%rho_r(i)%pw,error=error)
00313               END DO
00314               DO i=1,SIZE(epr_env%bind_set(idir,i_B)%rho%rho_g)
00315                  CALL pw_release(epr_env%bind_set(idir,i_B)%rho%rho_g(i)%pw,error=error)
00316               END DO
00317               DEALLOCATE(epr_env%bind_set(idir,i_B)%rho%rho_r,STAT=istat)
00318               CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00319               DEALLOCATE(epr_env%bind_set(idir,i_B)%rho%rho_g,STAT=istat)
00320               CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00321               DEALLOCATE(epr_env%bind_set(idir,i_B)%rho,STAT=istat)
00322               CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00323            ENDDO
00324            ENDDO
00325           DEALLOCATE(epr_env%bind_set,STAT=istat)
00326           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00327         END IF
00328         ! bind_atom_set
00329         IF(ASSOCIATED(epr_env%bind_atom_set)) THEN
00330           DEALLOCATE(epr_env%bind_atom_set,STAT=istat)
00331           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00332         END IF
00333         ! g_total
00334         IF(ASSOCIATED(epr_env%g_total)) THEN
00335           DEALLOCATE(epr_env%g_total,STAT=istat)
00336           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00337         END IF
00338         ! g_so
00339         IF(ASSOCIATED(epr_env%g_so)) THEN
00340           DEALLOCATE(epr_env%g_so,STAT=istat)
00341           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00342         END IF
00343         ! g_soo
00344         IF(ASSOCIATED(epr_env%g_soo)) THEN
00345           DEALLOCATE(epr_env%g_soo,STAT=istat)
00346           CPPostcondition(istat==0,cp_failure_level,routineP,error,failure)
00347         END IF
00348       END IF  ! ref count
00349     END IF ! failure
00350 
00351   END SUBROUTINE epr_env_cleanup
00352 
00353 END MODULE qs_linres_epr_utils