|
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 ! ***************************************************************************** 00011 MODULE paw_proj_set_types 00012 00013 USE basis_set_types, ONLY: get_gto_basis_set,& 00014 gto_basis_set_type 00015 USE cp_control_types, ONLY: qs_control_type 00016 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00017 cp_print_key_unit_nr 00018 USE input_section_types, ONLY: section_vals_type 00019 USE kinds, ONLY: default_string_length,& 00020 dp 00021 USE mathconstants, ONLY: dfac,& 00022 pi,& 00023 rootpi 00024 USE mathlib, ONLY: invert_matrix 00025 USE memory_utilities, ONLY: reallocate 00026 USE orbital_pointers, ONLY: indco,& 00027 indso,& 00028 nco,& 00029 ncoset,& 00030 nso,& 00031 nsoset 00032 USE orbital_transformation_matrices, ONLY: orbtramat 00033 USE qs_util, ONLY: exp_radius,& 00034 gauss_exponent 00035 USE termination, ONLY: stop_program 00036 #include "cp_common_uses.h" 00037 00038 IMPLICIT NONE 00039 00040 PRIVATE 00041 00042 ! *** Global parameters (only in this module) 00043 00044 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'paw_proj_set_types' 00045 00046 INTEGER, PARAMETER :: max_name_length = 60 00047 00048 ! *** Define the projectors types *** 00049 00050 ! ***************************************************************************** 00051 TYPE paw_proj_set_type 00052 INTEGER :: maxl,ncgauprj,nsgauprj 00053 INTEGER, DIMENSION(:), POINTER :: nprj ! 0:maxl 00054 INTEGER, DIMENSION(:), POINTER :: lx,ly,lz ! ncgauprj 00055 INTEGER, DIMENSION(:), POINTER :: ll, m ! nsgauprj 00056 INTEGER, DIMENSION(:), POINTER :: first_prj,last_prj ! 0:maxl 00057 INTEGER, DIMENSION(:), POINTER :: first_prjs ! 0:maxl 00058 REAL(dp) :: rcprj 00059 REAL(dp), DIMENSION(:), POINTER :: zisomin,zprjisomin 00060 REAL(dp), DIMENSION(:,:), POINTER :: zetprj ! maxnprj,0:maxl 00061 REAL(dp), DIMENSION(:,:), POINTER :: rzetprj ! maxnprj,0:maxl 00062 REAL(dp), DIMENSION(:,:), POINTER :: cprj ! ncgauprj, maxco*nset 00063 REAL(dp), DIMENSION(:,:), POINTER :: cprj_s ! nsgauprj, maxso*nset 00064 REAL(dp), DIMENSION(:,:), POINTER :: csprj ! ncgauprj, np_so 00065 REAL(dp), DIMENSION(:,:), POINTER :: local_oce_cphi_h,local_oce_cphi_s ! maxco,ncgf 00066 REAL(dp), DIMENSION(:,:), POINTER :: local_oce_sphi_h,local_oce_sphi_s ! maxco,nsgf 00067 REAL(dp), DIMENSION(:,:), POINTER :: sphi_h, sphi_s 00068 REAL(dp), DIMENSION(:,:,:,:), POINTER :: gccprj ! maxnprj,maxnpgf,0:maxl,nset 00069 LOGICAL, DIMENSION(:,:), POINTER :: isoprj ! maxnprj,0:maxl 00070 INTEGER :: nsatbas 00071 INTEGER :: nsotot 00072 INTEGER, DIMENSION(:), POINTER :: o2nindex ! maxso*nset 00073 INTEGER, DIMENSION(:), POINTER :: n2oindex ! maxso*nset 00074 00075 END TYPE paw_proj_set_type 00076 00077 ! *** Public subroutines *** 00078 00079 PUBLIC :: allocate_paw_proj_set,& 00080 deallocate_paw_proj_set,& 00081 get_paw_proj_set,& 00082 projectors,& 00083 set_paw_proj_set 00084 00085 ! *** Public data types *** 00086 00087 PUBLIC :: paw_proj_set_type 00088 00089 CONTAINS 00090 00091 ! ***************************************************************************** 00095 SUBROUTINE allocate_paw_proj_set(paw_proj_set,error) 00096 00097 TYPE(paw_proj_set_type), POINTER :: paw_proj_set 00098 TYPE(cp_error_type), INTENT(inout) :: error 00099 00100 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_paw_proj_set', 00101 routineP = moduleN//':'//routineN 00102 00103 INTEGER :: stat 00104 LOGICAL :: failure 00105 00106 failure = .FALSE. 00107 IF (ASSOCIATED(paw_proj_set)) CALL deallocate_paw_proj_set(paw_proj_set,error) 00108 00109 ALLOCATE (paw_proj_set,STAT=stat) 00110 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00111 00112 NULLIFY (paw_proj_set%nprj) 00113 NULLIFY (paw_proj_set%lx) 00114 NULLIFY (paw_proj_set%ly) 00115 NULLIFY (paw_proj_set%lz) 00116 NULLIFY (paw_proj_set%ll) 00117 NULLIFY (paw_proj_set%m) 00118 NULLIFY (paw_proj_set%first_prj) 00119 NULLIFY (paw_proj_set%last_prj) 00120 NULLIFY (paw_proj_set%first_prjs) 00121 00122 NULLIFY (paw_proj_set%zisomin) 00123 NULLIFY (paw_proj_set%zprjisomin) 00124 NULLIFY (paw_proj_set%zetprj) 00125 NULLIFY (paw_proj_set%cprj) 00126 NULLIFY (paw_proj_set%cprj_s) 00127 NULLIFY (paw_proj_set%csprj) 00128 NULLIFY (paw_proj_set%local_oce_cphi_h) 00129 NULLIFY (paw_proj_set%local_oce_cphi_s) 00130 NULLIFY (paw_proj_set%local_oce_sphi_h) 00131 NULLIFY (paw_proj_set%local_oce_sphi_s) 00132 NULLIFY (paw_proj_set%sphi_h) 00133 NULLIFY (paw_proj_set%sphi_s) 00134 NULLIFY (paw_proj_set%gccprj) 00135 NULLIFY (paw_proj_set%rzetprj) 00136 00137 NULLIFY (paw_proj_set%isoprj) 00138 00139 NULLIFY (paw_proj_set%o2nindex) 00140 NULLIFY (paw_proj_set%n2oindex) 00141 00142 END SUBROUTINE allocate_paw_proj_set 00143 00144 ! ***************************************************************************** 00148 SUBROUTINE deallocate_paw_proj_set(paw_proj_set,error) 00149 TYPE(paw_proj_set_type), POINTER :: paw_proj_set 00150 TYPE(cp_error_type), INTENT(inout) :: error 00151 00152 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_paw_proj_set', 00153 routineP = moduleN//':'//routineN 00154 00155 INTEGER :: stat 00156 LOGICAL :: failure 00157 00158 failure = .FALSE. 00159 00160 DEALLOCATE (paw_proj_set%zisomin,STAT=stat) 00161 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00162 DEALLOCATE (paw_proj_set%zprjisomin,STAT=stat) 00163 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00164 DEALLOCATE (paw_proj_set%nprj,STAT=stat) 00165 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00166 DEALLOCATE (paw_proj_set%lx,STAT=stat) 00167 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00168 DEALLOCATE (paw_proj_set%ly,STAT=stat) 00169 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00170 DEALLOCATE (paw_proj_set%lz,STAT=stat) 00171 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00172 DEALLOCATE (paw_proj_set%ll,STAT=stat) 00173 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00174 DEALLOCATE (paw_proj_set%m,STAT=stat) 00175 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00176 DEALLOCATE (paw_proj_set%first_prj,STAT=stat) 00177 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00178 DEALLOCATE (paw_proj_set%last_prj,STAT=stat) 00179 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00180 DEALLOCATE (paw_proj_set%first_prjs,STAT=stat) 00181 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00182 DEALLOCATE (paw_proj_set%zetprj,STAT=stat) 00183 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00184 DEALLOCATE (paw_proj_set%cprj,STAT=stat) 00185 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00186 DEALLOCATE (paw_proj_set%cprj_s,STAT=stat) 00187 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00188 DEALLOCATE (paw_proj_set%csprj,STAT=stat) 00189 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00190 DEALLOCATE (paw_proj_set%local_oce_cphi_h,STAT=stat) 00191 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00192 DEALLOCATE (paw_proj_set%local_oce_cphi_s,STAT=stat) 00193 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00194 DEALLOCATE (paw_proj_set%local_oce_sphi_h,STAT=stat) 00195 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00196 DEALLOCATE (paw_proj_set%local_oce_sphi_s,STAT=stat) 00197 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00198 DEALLOCATE (paw_proj_set%sphi_h,STAT=stat) 00199 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00200 DEALLOCATE (paw_proj_set%sphi_s,STAT=stat) 00201 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00202 DEALLOCATE (paw_proj_set%gccprj,STAT=stat) 00203 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00204 DEALLOCATE (paw_proj_set%isoprj,STAT=stat) 00205 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00206 DEALLOCATE (paw_proj_set%rzetprj,STAT=stat) 00207 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00208 DEALLOCATE (paw_proj_set%o2nindex,STAT=stat) 00209 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00210 DEALLOCATE (paw_proj_set%n2oindex,STAT=stat) 00211 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00212 00213 DEALLOCATE (paw_proj_set,STAT=stat) 00214 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00215 END SUBROUTINE deallocate_paw_proj_set 00216 00217 ! ***************************************************************************** 00221 SUBROUTINE projectors(paw_proj,orb_basis,rc,qs_control,max_rad_local_type,& 00222 force_env_section,error) 00223 00224 TYPE(paw_proj_set_type), POINTER :: paw_proj 00225 TYPE(gto_basis_set_type), POINTER :: orb_basis 00226 REAL(dp) :: rc 00227 TYPE(qs_control_type), INTENT(IN) :: qs_control 00228 REAL(dp), INTENT(IN) :: max_rad_local_type 00229 TYPE(section_vals_type), POINTER :: force_env_section 00230 TYPE(cp_error_type), INTENT(inout) :: error 00231 00232 CHARACTER(len=*), PARAMETER :: routineN = 'projectors', 00233 routineP = moduleN//':'//routineN 00234 00235 REAL(dp) :: eps_fit, eps_iso, eps_orb, 00236 eps_svd, max_rad_local 00237 00238 eps_fit = qs_control%gapw_control%eps_fit 00239 eps_iso = qs_control%gapw_control%eps_iso 00240 eps_svd = qs_control%gapw_control%eps_svd 00241 max_rad_local = qs_control%gapw_control%max_rad_local 00242 IF(max_rad_local_type.LT.max_rad_local)THEN 00243 max_rad_local=max_rad_local_type 00244 END IF 00245 eps_orb = qs_control%eps_pgf_orb 00246 00247 CALL build_projector(paw_proj,orb_basis,eps_fit,eps_iso,eps_svd, & 00248 rc,eps_orb,max_rad_local,force_env_section,error) 00249 00250 END SUBROUTINE projectors 00251 00252 ! ***************************************************************************** 00256 SUBROUTINE build_projector(paw_proj,orb_basis,eps_fit,eps_iso,eps_svd, & 00257 rc,eps_orb,max_rad_local,force_env_section,error) 00258 00259 TYPE(paw_proj_set_type), POINTER :: paw_proj 00260 TYPE(gto_basis_set_type), POINTER :: orb_basis 00261 REAL(dp), INTENT(IN) :: eps_fit, eps_iso, eps_svd, 00262 rc, eps_orb, max_rad_local 00263 TYPE(section_vals_type), POINTER :: force_env_section 00264 TYPE(cp_error_type), INTENT(inout) :: error 00265 00266 CHARACTER(len=*), PARAMETER :: routineN = 'build_projector', 00267 routineP = moduleN//':'//routineN 00268 00269 CHARACTER(LEN=default_string_length) :: bsname 00270 INTEGER :: ic, icgf, icgfmax, icgfmin, ico, icomax, icomin, il, info, ip, 00271 ipgf, ipp, iprj, iprjfirst, iprjs, is, iset, isgf, isgfmax, isgfmin, 00272 ishell, iso, iso_pgf, iso_set, isomin, jp, k, kp, lshell, lwork, lx, 00273 ly, lz, maxco, maxl, maxnprj, maxpgf, maxso, mp, ms, n, ncgauprj, ncgf, 00274 nisop, notprj, np, npgfg, ns, nset, nsgauprj, nsgf, nsox, output_unit, 00275 stat 00276 INTEGER, ALLOCATABLE, DIMENSION(:) :: IWORK 00277 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell 00278 INTEGER, DIMENSION(:, :), POINTER :: first_cgf, first_sgf, l, 00279 last_cgf, last_sgf 00280 LOGICAL :: failure 00281 LOGICAL, ALLOCATABLE, DIMENSION(:) :: isoprj 00282 REAL(dp) :: expzet, my_error, prefac, 00283 radius, x, zetmin, zetval 00284 REAL(dp), DIMENSION(:), POINTER :: zet, zetp 00285 REAL(dp), DIMENSION(:, :), POINTER :: cphi, gcc, smat, sphi, work 00286 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S, Work_dgesdd 00287 REAL(KIND=dp), ALLOCATABLE, 00288 DIMENSION(:, :) :: U, VT 00289 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius 00290 REAL(KIND=dp), DIMENSION(:, :, :), 00291 POINTER :: gcca, set_radius2 00292 TYPE(cp_logger_type), POINTER :: logger 00293 00294 NULLIFY(logger) 00295 logger => cp_error_get_logger(error) 00296 00297 NULLIFY (first_cgf,first_sgf,last_cgf,last_sgf,gcc,l,set_radius,set_radius2) 00298 NULLIFY (cphi,sphi,lmax,lmin,npgf,nshell,zet,zetp,smat,work,gcca) 00299 00300 failure =.FALSE. 00301 CPPrecondition(ASSOCIATED(paw_proj),cp_failure_level,routineP,error,failure) 00302 CPPrecondition(ASSOCIATED(orb_basis),cp_failure_level,routineP,error,failure) 00303 00304 IF(.NOT. failure) THEN 00305 CALL get_gto_basis_set(gto_basis_set=orb_basis,name=bsname,& 00306 ncgf=ncgf,nset=nset,nsgf=nsgf,& 00307 lmax=lmax,lmin=lmin,npgf=npgf,& 00308 nshell=nshell,cphi=cphi,sphi=sphi,& 00309 first_cgf=first_cgf,first_sgf=first_sgf,& 00310 l=l,last_cgf=last_cgf,last_sgf=last_sgf,& 00311 maxco=maxco,maxso=maxso,maxl=maxl,maxpgf=maxpgf,& 00312 gcc=gcca) 00313 00314 paw_proj%maxl = maxl 00315 CPPrecondition(.NOT. ASSOCIATED(paw_proj%zisomin),cp_failure_level,routineP,error,failure) 00316 CPPrecondition(.NOT. ASSOCIATED(paw_proj%zprjisomin),cp_failure_level,routineP,error,failure) 00317 CPPrecondition(.NOT. ASSOCIATED(paw_proj%nprj),cp_failure_level,routineP,error,failure) 00318 00319 ALLOCATE(paw_proj%zisomin(0:maxl),STAT=stat) 00320 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00321 paw_proj%zisomin(0:maxl) = 0.0_dp 00322 ALLOCATE(paw_proj%zprjisomin(0:maxl),STAT=stat) 00323 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00324 paw_proj%zprjisomin(0:maxl) = 0.0_dp 00325 ALLOCATE(paw_proj%nprj(0:maxl),STAT=stat) 00326 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00327 paw_proj%nprj(0:maxl) = 0 00328 00329 output_unit = cp_print_key_unit_nr(logger,force_env_section,& 00330 "DFT%PRINT%GAPW%PROJECTORS",extension=".Log",error=error) 00331 00332 IF (output_unit>0) THEN 00333 WRITE (UNIT=output_unit,FMT="(/,T2,A)")& 00334 "Projectors for the basis functions of "//TRIM(bsname) 00335 END IF 00336 00337 ALLOCATE(set_radius(nset),STAT=stat) 00338 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00339 set_radius = 0.0_dp 00340 DO iset = 1,nset 00341 DO is = 1,nshell(iset) 00342 set_radius(iset) = MAX(set_radius(iset), & 00343 exp_radius(l(is,iset),orb_basis%zet(npgf(iset),iset),eps_orb,gcca(npgf(iset),is,iset))) 00344 END DO ! is 00345 END DO ! iset 00346 00347 ALLOCATE(set_radius2(maxpgf,0:maxl,nset),STAT=stat) 00348 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00349 set_radius2 = 0.0_dp 00350 DO iset = 1,nset 00351 DO lshell = lmin(iset),lmax(iset) 00352 DO ip = 1, npgf(iset) 00353 set_radius2(ip,lshell,iset) = & 00354 exp_radius(lshell,orb_basis%zet(ip,iset),eps_orb,1.0_dp) 00355 END DO 00356 END DO ! is 00357 END DO ! iset 00358 00359 maxnprj = 0 00360 DO lshell = 0,maxl ! lshell 00361 np = 0 00362 DO iset = 1,nset 00363 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset))THEN 00364 DO ip = 1,npgf(iset) 00365 IF(set_radius2(ip,lshell,iset)<max_rad_local) THEN 00366 np = np + 1 00367 END IF 00368 END DO 00369 ENDIF 00370 ENDDO 00371 maxnprj = MAX(maxnprj,np) 00372 paw_proj%nprj(lshell) = np 00373 ENDDO ! lshell 00374 00375 ! *** Allocate exponents and coefficients *** 00376 ALLOCATE(paw_proj%zetprj(maxnprj,0:maxl),STAT=stat) 00377 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00378 paw_proj%zetprj(1:maxnprj,0:maxl)=0.0_dp 00379 ALLOCATE(paw_proj%gccprj(maxnprj,maxpgf,0:maxl,nset),STAT=stat) 00380 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00381 paw_proj%gccprj=0.0_dp 00382 ALLOCATE(paw_proj%rzetprj(maxnprj,0:maxl),STAT=stat) 00383 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00384 paw_proj%rzetprj(1:maxnprj,0:maxl)=0.0_dp 00385 ALLOCATE(paw_proj%isoprj(maxnprj,0:maxl),STAT=stat) 00386 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00387 paw_proj%isoprj = .FALSE. 00388 00389 mp = 0 00390 DO lshell = 0,maxl ! lshell 00391 mp = MAX(mp,paw_proj%nprj(lshell)*nco(lshell)) 00392 END DO 00393 00394 NULLIFY(zet,zetp,gcc,smat,work) 00395 ! *** Generate the projetor basis for each ang. mom. q.n. *** 00396 DO lshell = 0,maxl ! lshell 00397 00398 np = paw_proj%nprj(lshell) 00399 00400 ALLOCATE(isoprj(np),STAT=stat) 00401 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00402 isoprj = .FALSE. 00403 00404 ALLOCATE(zet(np),zetp(np),gcc(np,np),smat(np,np),work(np,np),STAT=stat) 00405 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00406 00407 zet(:) = 0.0_dp 00408 zetp(:) = 0.0_dp 00409 gcc(:,:) = 0.0_dp 00410 smat(:,:) = 0.0_dp 00411 work(:,:) = 0.0_dp 00412 00413 npgfg = 0 00414 notprj = 0 00415 ! *** Collect all the exponent which contribute to lshell *** 00416 DO iset = 1,nset ! iset 00417 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset))THEN 00418 DO ip = 1,npgf(iset) 00419 IF(set_radius2(ip,lshell,iset)<max_rad_local) THEN 00420 npgfg = npgfg + 1 00421 zet(npgfg) = orb_basis%zet(ip,iset) 00422 ELSE 00423 notprj = notprj + 1 00424 END IF 00425 END DO 00426 ENDIF 00427 ENDDO ! iset 00428 00429 ! *** Smallest exp. due to eps_iso: concerned as an isolated projector *** 00430 paw_proj%zisomin(lshell) = gauss_exponent(lshell,rc,eps_iso,1.0_dp) 00431 00432 nisop = 0 00433 DO ip = 1, np 00434 ! Check for equal exponents 00435 DO ipp =1,ip-1 00436 IF(zet(ip)==zet(ipp)) THEN 00437 CALL stop_program(routineP,moduleN,__LINE__,& 00438 "Linear dependency in the construction of the GAPW projectors:"//& 00439 " different sets of the BASIS SET contain identical exponents") 00440 END IF 00441 END DO 00442 00443 IF(zet(ip) >= paw_proj%zisomin(lshell)) THEN 00444 isoprj(ip) = .TRUE. 00445 nisop = nisop + 1 00446 ELSE 00447 isoprj(ip) = .FALSE. 00448 END IF 00449 END DO 00450 00451 ! *** Smallest exp. due to eps_fit: where to start geometric progression *** 00452 zetmin = gauss_exponent(lshell,rc,eps_fit,1.0_dp) 00453 00454 ! *** Generate the projectors by the geometric progression *** 00455 00456 IF(np-nisop-1 > 2) THEN 00457 x = (80.0_dp/zetmin)**(1.0_dp/REAL(np-nisop-1,dp)) 00458 ELSE 00459 x = 2.0_dp 00460 END IF 00461 IF(x > 2.0_dp) x = 2.0_dp 00462 00463 zetval = zetmin 00464 DO ip = np,1, -1 00465 IF(.NOT. isoprj(ip)) THEN 00466 zetp(ip) = zetval 00467 zetval = x*zetval 00468 ENDIF 00469 ENDDO 00470 00471 paw_proj%zprjisomin(lshell) = zetval 00472 00473 nisop = 0 00474 DO ip = np,1, -1 00475 IF(isoprj(ip)) THEN 00476 zetp(ip) = zetval 00477 zetval = x*zetval 00478 nisop = nisop + 1 00479 ENDIF 00480 ENDDO 00481 00482 ! *** Build the overlap matrix: <projector|primitive> *** 00483 prefac = 0.5_dp**(lshell + 2)*rootpi*dfac(2*lshell + 1) 00484 expzet = REAL(lshell,dp) + 1.5_dp 00485 00486 DO ip = 1,np 00487 IF (isoprj(ip)) THEN 00488 smat(ip,ip) = prefac/(zetp(ip) + zet(ip))**expzet 00489 ELSE 00490 DO jp = 1,np 00491 IF (.NOT. isoprj(jp)) THEN 00492 smat(ip,jp) = prefac/(zetp(ip) + zet(jp))**expzet 00493 ENDIF 00494 ENDDO 00495 ENDIF 00496 ENDDO 00497 00498 ! FIXME ... catch the case where np == 0 00499 CPPrecondition(np > 0,cp_failure_level,routineP,error,failure) 00500 00501 ! *** Compute inverse of the transpose *** 00502 IF (eps_svd.EQ.0.0_dp) THEN 00503 CALL invert_matrix(smat,gcc,my_error,"T",error=error) 00504 ELSE 00505 work=TRANSPOSE(smat) 00506 ! workspace query 00507 ALLOCATE(iwork(8*np),S(np),U(np,np),VT(np,np),work_dgesdd(1)) 00508 lwork=-1 00509 CALL DGESDD('S',np,np,work,np,S,U,np,vt,np,work_dgesdd,lwork,iwork,info) 00510 lwork=work_dgesdd(1) 00511 DEALLOCATE(work_dgesdd) ; ALLOCATE(work_dgesdd(lwork)) 00512 CALL DGESDD('S',np,np,work,np,S,U,np,vt,np,work_dgesdd,lwork,iwork,info) 00513 ! construct the inverse 00514 DO k=1,np 00515 ! invert SV 00516 IF (S(k)<eps_svd) THEN 00517 S(k)=0.0_dp 00518 ELSE 00519 S(k)=1.0_dp/S(k) 00520 ENDIF 00521 VT(k,:)=VT(k,:)*S(k) 00522 ENDDO 00523 CALL DGEMM('T','T',np,np,np,1.0_dp,VT,np,U,np,0.0_dp,gcc,np) 00524 DEALLOCATE(iwork,S,U,VT,work_dgesdd) 00525 ENDIF 00526 00527 ! *** Set the coefficient of the isolated projectors to 0 *** 00528 DO ip = 1, np 00529 IF(isoprj(ip)) THEN 00530 gcc(ip,ip) = 0.0_dp 00531 ENDIF 00532 ENDDO 00533 00534 ! *** Transfer data from local to global variables *** 00535 00536 paw_proj%zetprj(1:np,lshell) = zetp(1:np) 00537 paw_proj%isoprj(1:np,lshell) = isoprj(1:np) 00538 00539 npgfg = 0 00540 notprj = 0 00541 DO iset = 1,nset ! iset 00542 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset))THEN 00543 DO ip = 1,npgf(iset) 00544 IF(set_radius2(ip,lshell,iset)<max_rad_local) THEN 00545 npgfg = npgfg + 1 00546 paw_proj%gccprj(1:np,ip,lshell,iset) = gcc(1:np,npgfg) 00547 ELSE 00548 notprj = notprj + 1 00549 paw_proj%gccprj(1:np,ip,lshell,iset) = 0.0_dp 00550 END IF 00551 END DO 00552 ENDIF 00553 ENDDO ! iset 00554 00555 ! *** Print exponents and coefficients of the projectors *** 00556 IF (output_unit>0 ) THEN 00557 WRITE (UNIT=output_unit,FMT="(/,/,T2,A,I2)")& 00558 "Built projector for angular momentum quantum number l= ", lshell 00559 WRITE (UNIT=output_unit,FMT="(T2,A,I2)")& 00560 "Number of isolated projectors = ", nisop 00561 DO iset = 1,nset ! iset 00562 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset))THEN 00563 WRITE (UNIT=output_unit,FMT="(/,T2,A,I5,/,/,T4,A9,(T13,4f15.6))")& 00564 "Set ",iset,"exp prj: ", & 00565 (paw_proj%zetprj(ip,lshell),ip=1,np) 00566 DO jp = 1,npgf(iset) 00567 WRITE (UNIT=output_unit,FMT="(/,T4,A9,F15.6,/,T4,A9,(t13,4E15.6))") & 00568 "exp gto: ",orb_basis%zet(jp,iset),& 00569 "coeff.: ",(paw_proj%gccprj(ip,jp,& 00570 lshell,iset),ip=1,np) 00571 ENDDO 00572 ENDIF 00573 ENDDO ! iset 00574 END IF 00575 00576 ! *** Release the working storage for the current value lshell *** 00577 DEALLOCATE(isoprj,STAT=stat) 00578 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00579 DEALLOCATE(gcc,zet,zetp,smat,work,STAT=stat) 00580 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00581 00582 ENDDO ! lshell 00583 CALL cp_print_key_finished_output(output_unit,logger,force_env_section,& 00584 "DFT%PRINT%GAPW%PROJECTORS",error=error) 00585 00586 ! *** Release the working storage for the current value lshell *** 00587 DEALLOCATE(set_radius,STAT=stat) 00588 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00589 DEALLOCATE(set_radius2,STAT=stat) 00590 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00591 00592 ! *** Count primitives basis functions for the projectors 00593 paw_proj%ncgauprj = 0 00594 paw_proj%nsgauprj = 0 00595 DO lshell = 0,maxl 00596 paw_proj%ncgauprj = paw_proj%ncgauprj +nco(lshell)*paw_proj%nprj(lshell) 00597 paw_proj%nsgauprj = paw_proj%nsgauprj +nso(lshell)*paw_proj%nprj(lshell) 00598 ENDDO 00599 00600 ncgauprj = paw_proj%ncgauprj 00601 nsgauprj = paw_proj%nsgauprj 00602 CALL reallocate(paw_proj%cprj,1,ncgauprj,1,maxso*nset) 00603 CALL reallocate(paw_proj%lx,1,ncgauprj) 00604 CALL reallocate(paw_proj%ly,1,ncgauprj) 00605 CALL reallocate(paw_proj%lz,1,ncgauprj) 00606 CALL reallocate(paw_proj%first_prj,0,maxl) 00607 CALL reallocate(paw_proj%last_prj, 0,maxl) 00608 CALL reallocate(paw_proj%cprj_s,1,nsgauprj,1,maxso*nset) 00609 CALL reallocate(paw_proj%ll,1,nsgauprj) 00610 CALL reallocate(paw_proj%m,1,nsgauprj) 00611 CALL reallocate(paw_proj%first_prjs,0,maxl) 00612 00613 ncgauprj = 0 00614 nsgauprj = 0 00615 DO lshell = 0,maxl 00616 np = paw_proj%nprj(lshell) 00617 paw_proj%first_prj(lshell) = ncgauprj + 1 00618 paw_proj%first_prjs(lshell) = nsgauprj + 1 00619 paw_proj%last_prj(lshell) = ncgauprj + nco(lshell)*np 00620 DO ip = 1,np 00621 DO ico = ncoset(lshell-1) + 1, ncoset(lshell) 00622 ncgauprj = ncgauprj + 1 00623 paw_proj%lx(ncgauprj) = indco(1,ico) 00624 paw_proj%ly(ncgauprj) = indco(2,ico) 00625 paw_proj%lz(ncgauprj) = indco(3,ico) 00626 ENDDO ! ico 00627 DO iso = nsoset(lshell - 1) + 1 , nsoset(lshell) 00628 nsgauprj = nsgauprj + 1 00629 paw_proj%ll(nsgauprj) = indso(1,iso) 00630 paw_proj%m(nsgauprj) = indso(2,iso) 00631 ENDDO 00632 ENDDO ! ip 00633 ENDDO ! lshell 00634 00635 ms = 0 00636 DO iset=1,orb_basis%nset 00637 ns = nsoset(lmax(iset)) 00638 DO lshell = lmin(iset),lmax(iset) 00639 iprjfirst = paw_proj%first_prjs(lshell) 00640 np = paw_proj%nprj(lshell) 00641 DO ipgf = 1,npgf(iset) 00642 DO ip = 1, np 00643 DO il = 1, nso(lshell) 00644 iprjs = iprjfirst - 1 + il + (ip-1)* nso(lshell) 00645 iso = nsoset(lshell-1)+1+(lshell+paw_proj%m(iprjs)) 00646 00647 iso = iso + (ipgf-1)*ns + ms 00648 paw_proj%cprj_s(iprjs,iso) = & 00649 paw_proj%gccprj(ip,ipgf,lshell,iset) 00650 ENDDO ! iprjs 00651 ENDDO ! ip 00652 ENDDO ! ipgf 00653 ENDDO ! lshell 00654 ms = ms + maxso 00655 ENDDO ! iset 00656 00657 ms = 0 00658 DO iset=1,orb_basis%nset 00659 ns = nsoset(lmax(iset)) 00660 DO ipgf = 1,npgf(iset) 00661 DO iso = nsoset(lmin(iset)-1)+1+ns*(ipgf-1)+ms,ns*ipgf+ms 00662 DO lshell = 0,maxl 00663 np = paw_proj%nprj(lshell) 00664 DO ip = 1, np 00665 jp = paw_proj%first_prj(lshell) + (ip-1)* nco(lshell)-1 00666 kp = paw_proj%first_prjs(lshell)+ (ip-1)* nso(lshell)-1 00667 DO il = 1,nco(lshell) 00668 lx = indco(1,il+ncoset(lshell-1)) 00669 ly = indco(2,il+ncoset(lshell-1)) 00670 lz = indco(3,il+ncoset(lshell-1)) 00671 iprj = jp + il 00672 DO ic = 1,nso(lshell) 00673 iprjs = kp + ic 00674 paw_proj%cprj(iprj,iso) = paw_proj%cprj(iprj,iso) +& 00675 orbtramat(lshell)%s2c(ic,il)*paw_proj%cprj_s(iprjs,iso)*& 00676 SQRT(dfac(2*lshell+1)/(4.0_dp*pi)/& 00677 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1)) 00678 ENDDO ! ic 00679 ENDDO ! il 00680 ENDDO ! ip 00681 ENDDO ! lshell 00682 ENDDO ! iso 00683 ENDDO ! ipgf 00684 ms = ms + maxso 00685 ENDDO ! iset 00686 00687 ! local coefficients for the one center expansions : oce 00688 ! the coefficients are calculated for the full and soft expansions 00689 ALLOCATE(paw_proj%local_oce_cphi_h(maxco,ncgf),STAT=stat) 00690 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00691 paw_proj%local_oce_cphi_h = 0.0_dp 00692 ALLOCATE(paw_proj%local_oce_sphi_h(maxco,nsgf),STAT=stat) 00693 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00694 paw_proj%local_oce_sphi_h = 0.0_dp 00695 ALLOCATE(paw_proj%sphi_h(maxco,nsgf),STAT=stat) 00696 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00697 paw_proj%sphi_h = 0.0_dp 00698 00699 ALLOCATE(paw_proj%local_oce_cphi_s(maxco,ncgf),STAT=stat) 00700 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00701 paw_proj%local_oce_cphi_s = 0.0_dp 00702 ALLOCATE(paw_proj%local_oce_sphi_s(maxco,nsgf),STAT=stat) 00703 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00704 paw_proj%local_oce_sphi_s = 0.0_dp 00705 ALLOCATE(paw_proj%sphi_s(maxco,nsgf),STAT=stat) 00706 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00707 paw_proj%sphi_s = 0.0_dp 00708 00709 ! *** Cartesian *** 00710 DO iset = 1,nset 00711 n = ncoset(lmax(iset)) 00712 DO ipgf = 1,npgf(iset) 00713 DO ishell = 1,nshell(iset) 00714 lshell = l(ishell,iset) 00715 icomin = ncoset(lshell-1) + 1 + n*(ipgf - 1) 00716 icomax = ncoset(lshell) + n*(ipgf - 1) 00717 icgfmin = first_cgf(ishell,iset) 00718 icgfmax = last_cgf(ishell,iset) 00719 radius = exp_radius(lshell,orb_basis%zet(ipgf,iset),& 00720 eps_fit,1.0_dp) 00721 DO icgf = icgfmin,icgfmax 00722 paw_proj%local_oce_cphi_h(icomin:icomax,icgf) = & 00723 cphi(icomin:icomax,icgf) 00724 IF (radius < rc ) THEN 00725 paw_proj%local_oce_cphi_s(icomin:icomax,icgf)=0.0_dp 00726 ELSE 00727 paw_proj%local_oce_cphi_s(icomin:icomax,icgf) = & 00728 cphi(icomin:icomax,icgf) 00729 ENDIF 00730 END DO 00731 ENDDO ! ishell 00732 ENDDO ! ipgf 00733 ENDDO ! iset 00734 00735 DO iset = 1,nset 00736 n = ncoset(lmax(iset)) 00737 DO ipgf = 1,npgf(iset) 00738 DO ishell = 1,nshell(iset) 00739 lshell = l(ishell,iset) 00740 icomin = ncoset(lshell-1) + 1 + n*(ipgf - 1) 00741 icomax = ncoset(lshell) + n*(ipgf - 1) 00742 isgfmin = first_sgf(ishell,iset) 00743 isgfmax = last_sgf(ishell,iset) 00744 radius = exp_radius(lshell,orb_basis%zet(ipgf,iset),& 00745 eps_fit,1.0_dp) 00746 DO isgf = isgfmin,isgfmax 00747 paw_proj%sphi_h(icomin:icomax,isgf) = & 00748 sphi(icomin:icomax,isgf) 00749 IF (radius < rc ) THEN 00750 paw_proj%sphi_s(icomin:icomax,isgf)=0.0_dp 00751 ELSE 00752 paw_proj%sphi_s(icomin:icomax,isgf) = & 00753 sphi(icomin:icomax,isgf) 00754 ENDIF 00755 END DO 00756 ENDDO ! ishell 00757 ENDDO ! ipgf 00758 ENDDO ! iset 00759 00760 DO iset = 1,nset 00761 n = ncoset(lmax(iset)) 00762 ns = nsoset(lmax(iset)) 00763 DO ipgf = 1,npgf(iset) 00764 DO ishell = 1,nshell(iset) 00765 lshell = l(ishell,iset) 00766 icomin = ncoset(lshell-1) + 1 + n*(ipgf - 1) 00767 icomax = ncoset(lshell) + n*(ipgf - 1) 00768 isgfmin = first_sgf(ishell,iset) 00769 isgfmax = last_sgf(ishell,iset) 00770 isomin = nsoset(lshell-1) + 1 +ns*(ipgf - 1) 00771 DO is = 1,nso(lshell) 00772 iso = isomin + is - 1 00773 DO ic = 1,nco(lshell) 00774 ico = icomin + ic -1 00775 lx = indco(1,ic+ncoset(lshell-1)) 00776 ly = indco(2,ic+ncoset(lshell-1)) 00777 lz = indco(3,ic+ncoset(lshell-1)) 00778 DO isgf = isgfmin,isgfmax 00779 paw_proj%local_oce_sphi_h(iso,isgf) = & 00780 paw_proj%local_oce_sphi_h(iso,isgf) + & 00781 orbtramat(lshell)%s2c(is,ic)*& 00782 paw_proj%sphi_h(ico,isgf)*SQRT((4.0_dp*pi*& 00783 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1))/dfac(2*lshell+1)) 00784 paw_proj%local_oce_sphi_s(iso,isgf) = & 00785 paw_proj%local_oce_sphi_s(iso,isgf) + & 00786 orbtramat(lshell)%s2c(is,ic)*& 00787 paw_proj%sphi_s(ico,isgf)*SQRT((4.0_dp*pi*& 00788 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1))/dfac(2*lshell+1)) 00789 END DO ! isgf 00790 END DO ! ic 00791 END DO ! is 00792 END DO ! ishell 00793 END DO ! ipgf 00794 END DO ! iset 00795 00796 !Index transformation OLD-NEW 00797 ALLOCATE(paw_proj%o2nindex(maxso*nset),STAT=stat) 00798 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00799 ALLOCATE(paw_proj%n2oindex(maxso*nset),STAT=stat) 00800 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00801 paw_proj%o2nindex = 0 00802 paw_proj%n2oindex = 0 00803 ico = 1 00804 DO iset=1,nset 00805 iso_set = (iset - 1)*maxso + 1 00806 nsox = nsoset(lmax(iset)) 00807 DO ipgf=1,npgf(iset) 00808 iso_pgf = iso_set + (ipgf - 1)*nsox 00809 iso = iso_pgf + nsoset(lmin(iset)-1) 00810 DO lx=lmin(iset),lmax(iset) 00811 DO k=1,nso(lx) 00812 paw_proj%n2oindex(ico) = iso 00813 paw_proj%o2nindex(iso) = ico 00814 iso = iso + 1 00815 ico = ico + 1 00816 END DO 00817 END DO 00818 END DO 00819 END DO 00820 mp = ico-1 00821 paw_proj%nsatbas = mp 00822 paw_proj%nsotot = nset*maxso 00823 ALLOCATE(paw_proj%csprj(nsgauprj,mp),STAT=stat) 00824 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00825 paw_proj%csprj = 0._dp 00826 DO k=1,mp 00827 ico = paw_proj%n2oindex(k) 00828 paw_proj%csprj(:,k) = paw_proj%cprj_s(:,ico) 00829 END DO 00830 00831 END IF ! failure 00832 END SUBROUTINE build_projector 00833 00834 ! ***************************************************************************** 00838 SUBROUTINE get_paw_proj_set(paw_proj_set,cprj,cprj_s,csprj,& 00839 first_prj,first_prjs,last_prj, & 00840 local_oce_sphi_h,local_oce_sphi_s,& 00841 maxl,ncgauprj,nsgauprj,nsatbas,nsotot,nprj,& 00842 o2nindex,n2oindex,& 00843 rcprj,rzetprj,zisomin,zetprj) 00844 00845 TYPE(paw_proj_set_type), POINTER :: paw_proj_set 00846 REAL(dp), DIMENSION(:, :), OPTIONAL, 00847 POINTER :: cprj, cprj_s, csprj 00848 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: first_prj, first_prjs, 00849 last_prj 00850 REAL(dp), DIMENSION(:, :), OPTIONAL, 00851 POINTER :: local_oce_sphi_h, 00852 local_oce_sphi_s 00853 INTEGER, INTENT(OUT), OPTIONAL :: maxl, ncgauprj, nsgauprj, 00854 nsatbas, nsotot 00855 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nprj, o2nindex, n2oindex 00856 REAL(dp), INTENT(OUT), OPTIONAL :: rcprj 00857 REAL(dp), DIMENSION(:, :), OPTIONAL, 00858 POINTER :: rzetprj 00859 REAL(dp), DIMENSION(:), OPTIONAL, 00860 POINTER :: zisomin 00861 REAL(dp), DIMENSION(:, :), OPTIONAL, 00862 POINTER :: zetprj 00863 00864 CHARACTER(len=*), PARAMETER :: routineN = 'get_paw_proj_set', 00865 routineP = moduleN//':'//routineN 00866 00867 IF (ASSOCIATED(paw_proj_set)) THEN 00868 IF (PRESENT(cprj)) cprj => paw_proj_set%cprj 00869 IF (PRESENT(cprj_s)) cprj_s => paw_proj_set%cprj_s 00870 IF (PRESENT(csprj)) csprj => paw_proj_set%csprj 00871 IF (PRESENT(local_oce_sphi_h)) local_oce_sphi_h => paw_proj_set%local_oce_sphi_h 00872 IF (PRESENT(local_oce_sphi_s)) local_oce_sphi_s => paw_proj_set%local_oce_sphi_s 00873 IF (PRESENT(first_prj)) first_prj => paw_proj_set%first_prj 00874 IF (PRESENT(last_prj)) last_prj => paw_proj_set%last_prj 00875 IF (PRESENT(first_prjs)) first_prjs => paw_proj_set%first_prjs 00876 IF (PRESENT(maxl)) maxl = paw_proj_set%maxl 00877 IF (PRESENT(ncgauprj )) ncgauprj = paw_proj_set%ncgauprj 00878 IF (PRESENT(nsgauprj )) nsgauprj = paw_proj_set%nsgauprj 00879 IF (PRESENT(nsatbas )) nsatbas = paw_proj_set%nsatbas 00880 IF (PRESENT(nsotot )) nsotot = paw_proj_set%nsotot 00881 IF (PRESENT(nprj)) nprj => paw_proj_set%nprj 00882 IF (PRESENT(rcprj)) rcprj = paw_proj_set%rcprj 00883 IF (PRESENT(rzetprj)) rzetprj => paw_proj_set%rzetprj 00884 IF (PRESENT(zisomin)) zisomin => paw_proj_set%zisomin 00885 IF (PRESENT(zetprj)) zetprj => paw_proj_set%zetprj 00886 IF (PRESENT(o2nindex)) o2nindex => paw_proj_set%o2nindex 00887 IF (PRESENT(n2oindex)) n2oindex => paw_proj_set%n2oindex 00888 00889 ELSE 00890 CALL stop_program(routineP,moduleN,__LINE__,& 00891 "The pointer paw_proj_set is not associated") 00892 ENDIF 00893 00894 END SUBROUTINE get_paw_proj_set 00895 00896 ! ***************************************************************************** 00900 SUBROUTINE set_paw_proj_set(paw_proj_set, rzetprj, rcprj) 00901 00902 TYPE(paw_proj_set_type), POINTER :: paw_proj_set 00903 REAL(dp), DIMENSION(:, :), OPTIONAL, 00904 POINTER :: rzetprj 00905 REAL(dp), INTENT(IN), OPTIONAL :: rcprj 00906 00907 CHARACTER(len=*), PARAMETER :: routineN = 'set_paw_proj_set', 00908 routineP = moduleN//':'//routineN 00909 00910 IF (ASSOCIATED(paw_proj_set)) THEN 00911 IF (PRESENT(rzetprj)) paw_proj_set%rzetprj(:,0:) = rzetprj(:,0:) 00912 IF (PRESENT(rcprj)) paw_proj_set%rcprj = rcprj 00913 ELSE 00914 CALL stop_program(routineP,moduleN,__LINE__,& 00915 "The pointer paw_proj_set is not associated") 00916 END IF 00917 00918 END SUBROUTINE set_paw_proj_set 00919 00920 END MODULE paw_proj_set_types 00921
1.7.3