|
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 ! ***************************************************************************** 00013 MODULE atom_operators 00014 USE ai_onecenter, ONLY: & 00015 sg_conf, sg_coulomb, sg_erf, sg_exchange, sg_gpot, sg_kinetic, & 00016 sg_kinnuc, sg_nuclear, sg_overlap, sg_proj_ol, sto_conf, sto_kinetic, & 00017 sto_nuclear, sto_overlap 00018 USE atom_types, ONLY: & 00019 CGTO_BASIS, GTH_PSEUDO, GTO_BASIS, NO_PSEUDO, NUM_BASIS, STO_BASIS, & 00020 atom_basis_type, atom_integrals, atom_potential_type, atom_state 00021 USE atom_utils, ONLY: coulomb_potential_numeric,& 00022 integrate_grid,& 00023 numpot_matrix,& 00024 slater_density,& 00025 wigner_slater_functional 00026 USE dkh_main, ONLY: dkh_atom_transformation 00027 USE erf_fn, ONLY: erfc 00028 USE f77_blas 00029 USE input_constants, ONLY: do_dkh0_atom,& 00030 do_dkh1_atom,& 00031 do_dkh2_atom,& 00032 do_dkh3_atom,& 00033 do_dkh4_atom,& 00034 do_dkh5_atom,& 00035 do_nonrel_atom,& 00036 do_zoramp_atom 00037 USE kinds, ONLY: dp 00038 USE lapack, ONLY: lapack_sgesv,& 00039 lapack_ssyev 00040 USE mathconstants, ONLY: gamma1,& 00041 sqrt2 00042 USE periodic_table, ONLY: ptable 00043 USE physcon, ONLY: a_fine 00044 USE qs_grid_atom, ONLY: grid_atom_type 00045 USE timings, ONLY: timeset,& 00046 timestop 00047 #include "cp_common_uses.h" 00048 00049 IMPLICIT NONE 00050 00051 PRIVATE 00052 00053 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_operators' 00054 00055 PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release 00056 PUBLIC :: atom_relint_setup, atom_relint_release 00057 00058 ! ***************************************************************************** 00059 00060 CONTAINS 00061 00062 ! ***************************************************************************** 00063 SUBROUTINE atom_int_setup(integrals,basis,potential,& 00064 eri_coulomb,eri_exchange,all_nu,error) 00065 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00066 TYPE(atom_basis_type), INTENT(INOUT) :: basis 00067 TYPE(atom_potential_type), INTENT(IN), 00068 OPTIONAL :: potential 00069 LOGICAL, INTENT(IN), OPTIONAL :: eri_coulomb, eri_exchange, 00070 all_nu 00071 TYPE(cp_error_type), INTENT(INOUT) :: error 00072 00073 CHARACTER(len=*), PARAMETER :: routineN = 'atom_int_setup', 00074 routineP = moduleN//':'//routineN 00075 00076 INTEGER :: handle, i, ierr, ii, info, 00077 ipiv(1000), l, l1, l2, ll, 00078 lwork, m, m1, m2, mm1, mm2, 00079 n, n1, n2, nn1, nn2, nu, nx 00080 LOGICAL :: failure 00081 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work 00082 REAL(KIND=dp), ALLOCATABLE, 00083 DIMENSION(:, :) :: omat 00084 REAL(KIND=dp), DIMENSION(:, :), POINTER :: eri 00085 00086 failure = .FALSE. 00087 00088 CALL timeset(routineN,handle) 00089 00090 IF ( integrals%status == 0 ) THEN 00091 n = MAXVAL(basis%nbas) 00092 integrals%n = basis%nbas 00093 00094 IF ( PRESENT(eri_coulomb) ) THEN 00095 integrals%eri_coulomb = eri_coulomb 00096 ELSE 00097 integrals%eri_coulomb = .FALSE. 00098 END IF 00099 IF ( PRESENT(eri_exchange) ) THEN 00100 integrals%eri_exchange = eri_exchange 00101 ELSE 00102 integrals%eri_exchange = .FALSE. 00103 END IF 00104 IF ( PRESENT(all_nu) ) THEN 00105 integrals%all_nu = all_nu 00106 ELSE 00107 integrals%all_nu = .FALSE. 00108 END IF 00109 00110 NULLIFY ( integrals%ovlp, integrals%kin, integrals%core, integrals%conf ) 00111 DO ll=1,SIZE(integrals%ceri) 00112 NULLIFY ( integrals%ceri(ll)%int, integrals%eeri(ll)%int ) 00113 END DO 00114 00115 ALLOCATE (integrals%ovlp(n,n,0:3),STAT=ierr) 00116 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00117 integrals%ovlp = 0._dp 00118 00119 ALLOCATE (integrals%kin(n,n,0:3),STAT=ierr) 00120 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00121 integrals%kin = 0._dp 00122 00123 integrals%status = 1 00124 00125 IF ( PRESENT(potential) ) THEN 00126 IF ( potential%confinement ) THEN 00127 ALLOCATE (integrals%conf(n,n,0:3),STAT=ierr) 00128 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00129 integrals%conf = 0._dp 00130 END IF 00131 END IF 00132 00133 SELECT CASE (basis%basis_type) 00134 CASE DEFAULT 00135 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00136 CASE (GTO_BASIS) 00137 DO l=0,3 00138 n = integrals%n(l) 00139 CALL sg_overlap ( integrals%ovlp(1:n,1:n,l), l, basis%am(1:n,l), basis%am(1:n,l) ) 00140 CALL sg_kinetic ( integrals%kin(1:n,1:n,l), l, basis%am(1:n,l), basis%am(1:n,l) ) 00141 IF ( PRESENT(potential) ) THEN 00142 IF ( potential%confinement ) THEN 00143 CALL sg_conf ( integrals%conf(1:n,1:n,l), potential%rcon, potential%ncon/2, & 00144 l, basis%am(1:n,l), basis%am(1:n,l) ) 00145 END IF 00146 END IF 00147 END DO 00148 IF ( integrals%eri_coulomb ) THEN 00149 ll = 0 00150 DO l1=0,3 00151 n1 = integrals%n(l1) 00152 nn1 = (n1*(n1+1))/2 00153 DO l2=0,l1 00154 n2 = integrals%n(l2) 00155 nn2 = (n2*(n2+1))/2 00156 IF ( integrals%all_nu ) THEN 00157 nx = MIN(2*l1,2*l2) 00158 ELSE 00159 nx = 0 00160 END IF 00161 DO nu = 0, nx, 2 00162 ll = ll + 1 00163 ALLOCATE (integrals%ceri(ll)%int(nn1,nn2),STAT=ierr) 00164 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00165 integrals%ceri(ll)%int = 0._dp 00166 eri => integrals%ceri(ll)%int 00167 CALL sg_coulomb ( eri, nu, basis%am(1:n1,l1), l1, basis%am(1:n2,l2), l2 ) 00168 END DO 00169 END DO 00170 END DO 00171 END IF 00172 IF ( integrals%eri_exchange ) THEN 00173 ll = 0 00174 DO l1=0,3 00175 n1 = integrals%n(l1) 00176 nn1 = (n1*(n1+1))/2 00177 DO l2=0,l1 00178 n2 = integrals%n(l2) 00179 nn2 = (n2*(n2+1))/2 00180 DO nu = ABS(l1-l2),l1+l2,2 00181 ll = ll + 1 00182 ALLOCATE (integrals%eeri(ll)%int(nn1,nn2),STAT=ierr) 00183 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00184 integrals%eeri(ll)%int = 0._dp 00185 eri => integrals%eeri(ll)%int 00186 CALL sg_exchange ( eri, nu, basis%am(1:n1,l1), l1, basis%am(1:n2,l2), l2 ) 00187 END DO 00188 END DO 00189 END DO 00190 END IF 00191 CASE (CGTO_BASIS) 00192 DO l=0,3 00193 n = integrals%n(l) 00194 m = basis%nprim(l) 00195 IF (n>0 .AND. m>0) THEN 00196 ALLOCATE (omat(m,m),STAT=ierr) 00197 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00198 00199 CALL sg_overlap ( omat(1:m,1:m), l, basis%am(1:m,l), basis%am(1:m,l) ) 00200 CALL contract2(integrals%ovlp(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00201 CALL sg_kinetic ( omat(1:m,1:m), l, basis%am(1:m,l), basis%am(1:m,l) ) 00202 CALL contract2(integrals%kin(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00203 IF ( PRESENT(potential) ) THEN 00204 IF ( potential%confinement ) THEN 00205 CALL sg_conf ( omat(1:m,1:m), potential%rcon, potential%ncon/2, & 00206 l, basis%am(1:m,l), basis%am(1:m,l) ) 00207 CALL contract2(integrals%conf(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00208 END IF 00209 END IF 00210 DEALLOCATE (omat,STAT=ierr) 00211 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00212 END IF 00213 END DO 00214 IF ( integrals%eri_coulomb ) THEN 00215 ll = 0 00216 DO l1=0,3 00217 n1 = integrals%n(l1) 00218 nn1 = (n1*(n1+1))/2 00219 m1 = basis%nprim(l1) 00220 mm1 = (m1*(m1+1))/2 00221 DO l2=0,l1 00222 n2 = integrals%n(l2) 00223 nn2 = (n2*(n2+1))/2 00224 m2 = basis%nprim(l2) 00225 mm2 = (m2*(m2+1))/2 00226 IF ( integrals%all_nu ) THEN 00227 nx = MIN(2*l1,2*l2) 00228 ELSE 00229 nx = 0 00230 END IF 00231 DO nu = 0, nx, 2 00232 ll = ll + 1 00233 ALLOCATE (integrals%ceri(ll)%int(nn1,nn2),STAT=ierr) 00234 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00235 integrals%ceri(ll)%int = 0._dp 00236 ALLOCATE (omat(mm1,mm2),STAT=ierr) 00237 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00238 00239 eri => integrals%ceri(ll)%int 00240 CALL sg_coulomb ( omat, nu, basis%am(1:m1,l1), l1, basis%am(1:m2,l2), l2 ) 00241 CALL contract4 ( eri, omat, basis%cm(1:m1,1:n1,l1), basis%cm(1:m2,1:n2,l2), error ) 00242 00243 DEALLOCATE (omat,STAT=ierr) 00244 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00245 END DO 00246 END DO 00247 END DO 00248 END IF 00249 IF ( integrals%eri_exchange ) THEN 00250 ll = 0 00251 DO l1=0,3 00252 n1 = integrals%n(l1) 00253 nn1 = (n1*(n1+1))/2 00254 m1 = basis%nprim(l1) 00255 mm1 = (m1*(m1+1))/2 00256 DO l2=0,l1 00257 n2 = integrals%n(l2) 00258 nn2 = (n2*(n2+1))/2 00259 m2 = basis%nprim(l2) 00260 mm2 = (m2*(m2+1))/2 00261 DO nu = ABS(l1-l2),l1+l2,2 00262 ll = ll + 1 00263 ALLOCATE (integrals%eeri(ll)%int(nn1,nn2),STAT=ierr) 00264 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00265 integrals%eeri(ll)%int = 0._dp 00266 ALLOCATE (omat(mm1,mm2),STAT=ierr) 00267 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00268 00269 eri => integrals%eeri(ll)%int 00270 CALL sg_exchange ( omat, nu, basis%am(1:m1,l1), l1, basis%am(1:m2,l2), l2 ) 00271 CALL contract4 ( eri, omat, basis%cm(1:m1,1:n1,l1), basis%cm(1:m2,1:n2,l2), error ) 00272 00273 DEALLOCATE (omat,STAT=ierr) 00274 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00275 END DO 00276 END DO 00277 END DO 00278 END IF 00279 CASE (STO_BASIS) 00280 DO l=0,3 00281 n = integrals%n(l) 00282 CALL sto_overlap ( integrals%ovlp(1:n,1:n,l), basis%ns(1:n,l), basis%as(1:n,l), & 00283 basis%ns(1:n,l), basis%as(1:n,l) ) 00284 CALL sto_kinetic ( integrals%kin(1:n,1:n,l), l, basis%ns(1:n,l), basis%as(1:n,l), & 00285 basis%ns(1:n,l), basis%as(1:n,l) ) 00286 IF ( PRESENT(potential) ) THEN 00287 IF ( potential%confinement ) THEN 00288 CALL sto_conf ( integrals%conf(1:n,1:n,l), potential%rcon, potential%ncon/2, & 00289 basis%ns(1:n,l), basis%as(1:n,l), basis%ns(1:n,l), basis%as(1:n,l) ) 00290 END IF 00291 END IF 00292 END DO 00293 CPAssert(.NOT.integrals%eri_coulomb,cp_failure_level,routineP,error,failure) 00294 CPAssert(.NOT.integrals%eri_exchange,cp_failure_level,routineP,error,failure) 00295 CASE (NUM_BASIS) 00296 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00297 END SELECT 00298 00299 ! setup transformation matrix to get an orthogonal basis, remove linear dependencies 00300 NULLIFY(integrals%utrans,integrals%uptrans) 00301 n = MAXVAL(basis%nbas) 00302 ALLOCATE (integrals%utrans(n,n,0:3),integrals%uptrans(n,n,0:3),STAT=ierr) 00303 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00304 integrals%utrans = 0._dp 00305 integrals%uptrans = 0._dp 00306 integrals%nne = integrals%n 00307 lwork = 10*n 00308 ALLOCATE(omat(n,n),w(n),work(lwork),STAT=ierr) 00309 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00310 DO l = 0, 3 00311 n = integrals%n(l) 00312 IF ( n > 0 ) THEN 00313 omat(1:n,1:n) = integrals%ovlp(1:n,1:n,l) 00314 CALL lapack_ssyev ( "V", "U", n, omat(1:n,1:n), n, w(1:n), work, lwork, info ) 00315 CPPostcondition(info==0,cp_failure_level,routineP,error,failure) 00316 ii = 0 00317 DO i=1,n 00318 IF (w(i) > basis%eps_eig) THEN 00319 ii = ii + 1 00320 integrals%utrans(1:n,ii,l) = omat(1:n,i)/SQRT(w(i)) 00321 END IF 00322 END DO 00323 integrals%nne(l) = ii 00324 omat(1:ii,1:ii)=MATMUL(TRANSPOSE(integrals%utrans(1:n,1:ii,l)),integrals%utrans(1:n,1:ii,l)) 00325 DO i=1,ii 00326 integrals%uptrans(i,i,l)=1._dp 00327 ENDDO 00328 CALL lapack_sgesv ( ii, ii, omat(1:ii,1:ii), ii, ipiv, integrals%uptrans(1:ii,1:ii,l), ii, info ) 00329 CPPostcondition(info==0,cp_failure_level,routineP,error,failure) 00330 END IF 00331 END DO 00332 DEALLOCATE(omat,w,work,STAT=ierr) 00333 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00334 00335 END IF 00336 00337 CALL timestop(handle) 00338 00339 END SUBROUTINE atom_int_setup 00340 ! ***************************************************************************** 00341 SUBROUTINE atom_ppint_setup(integrals,basis,potential,error) 00342 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00343 TYPE(atom_basis_type), INTENT(INOUT) :: basis 00344 TYPE(atom_potential_type), INTENT(IN) :: potential 00345 TYPE(cp_error_type), INTENT(INOUT) :: error 00346 00347 CHARACTER(len=*), PARAMETER :: routineN = 'atom_ppint_setup', 00348 routineP = moduleN//':'//routineN 00349 00350 INTEGER :: handle, i, ierr, ii, j, k, l, 00351 m, n 00352 LOGICAL :: failure 00353 REAL(KIND=dp) :: al, alpha, rc 00354 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, xmat 00355 REAL(KIND=dp), ALLOCATABLE, 00356 DIMENSION(:, :) :: omat, spmat 00357 REAL(KIND=dp), DIMENSION(:), POINTER :: rad 00358 00359 failure = .FALSE. 00360 00361 CALL timeset(routineN,handle) 00362 00363 IF ( integrals%ppstat == 0 ) THEN 00364 n = MAXVAL(basis%nbas) 00365 integrals%n = basis%nbas 00366 00367 NULLIFY ( integrals%core, integrals%hnl ) 00368 00369 ALLOCATE (integrals%hnl(n,n,0:3),STAT=ierr) 00370 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00371 integrals%hnl = 0._dp 00372 00373 ALLOCATE (integrals%core(n,n,0:3),STAT=ierr) 00374 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00375 integrals%core = 0._dp 00376 00377 ALLOCATE (integrals%clsd(n,n,0:3),STAT=ierr) 00378 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00379 integrals%clsd = 0._dp 00380 00381 integrals%ppstat = 1 00382 00383 SELECT CASE (basis%basis_type) 00384 CASE DEFAULT 00385 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00386 CASE (GTO_BASIS) 00387 00388 SELECT CASE (potential%ppot_type) 00389 CASE (NO_PSEUDO) 00390 DO l=0,3 00391 n = integrals%n(l) 00392 CALL sg_nuclear ( integrals%core(1:n,1:n,l), l, basis%am(1:n,l), basis%am(1:n,l) ) 00393 END DO 00394 CASE (GTH_PSEUDO) 00395 alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp) 00396 DO l=0,3 00397 n = integrals%n(l) 00398 ALLOCATE (omat(n,n),spmat(n,5),STAT=ierr) 00399 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00400 00401 omat = 0._dp 00402 CALL sg_erf ( omat(1:n,1:n), l, alpha, basis%am(1:n,l), basis%am(1:n,l) ) 00403 integrals%core(1:n,1:n,l) = -potential%gth_pot%zion*omat(1:n,1:n) 00404 DO i=1,potential%gth_pot%ncl 00405 omat = 0._dp 00406 CALL sg_gpot ( omat(1:n,1:n), i-1, potential%gth_pot%rc, l, basis%am(1:n,l), basis%am(1:n,l) ) 00407 integrals%core(1:n,1:n,l) = integrals%core(1:n,1:n,l) + & 00408 potential%gth_pot%cl(i)*omat(1:n,1:n) 00409 END DO 00410 IF (potential%gth_pot%lpotextended) THEN 00411 DO k=1,potential%gth_pot%nexp_lpot 00412 DO i=1,potential%gth_pot%nct_lpot(k) 00413 omat = 0._dp 00414 CALL sg_gpot ( omat(1:n,1:n), i-1, potential%gth_pot%alpha_lpot(k), l, & 00415 basis%am(1:n,l), basis%am(1:n,l) ) 00416 integrals%core(1:n,1:n,l) = integrals%core(1:n,1:n,l) + & 00417 potential%gth_pot%cval_lpot(i,k)*omat(1:n,1:n) 00418 END DO 00419 END DO 00420 END IF 00421 IF (potential%gth_pot%lsdpot) THEN 00422 DO k=1,potential%gth_pot%nexp_lsd 00423 DO i=1,potential%gth_pot%nct_lsd(k) 00424 omat = 0._dp 00425 CALL sg_gpot ( omat(1:n,1:n), i-1, potential%gth_pot%alpha_lsd(k), l, & 00426 basis%am(1:n,l), basis%am(1:n,l) ) 00427 integrals%clsd(1:n,1:n,l) = integrals%clsd(1:n,1:n,l) + & 00428 potential%gth_pot%cval_lsd(i,k)*omat(1:n,1:n) 00429 END DO 00430 END DO 00431 END IF 00432 00433 spmat = 0._dp 00434 m = potential%gth_pot%nl(l) 00435 DO i=1,m 00436 CALL sg_proj_ol ( spmat(1:n,i), l, basis%am(1:n,l), i-1, potential%gth_pot%rcnl(l) ) 00437 END DO 00438 integrals%hnl(1:n,1:n,l) = MATMUL(spmat(1:n,1:m),& 00439 MATMUL(potential%gth_pot%hnl(1:m,1:m,l),TRANSPOSE(spmat(1:n,1:m)))) 00440 00441 DEALLOCATE (omat,spmat,STAT=ierr) 00442 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00443 END DO 00444 CASE DEFAULT 00445 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00446 END SELECT 00447 00448 CASE (CGTO_BASIS) 00449 00450 SELECT CASE (potential%ppot_type) 00451 CASE (NO_PSEUDO) 00452 DO l=0,3 00453 n = integrals%n(l) 00454 m = basis%nprim(l) 00455 ALLOCATE (omat(m,m),STAT=ierr) 00456 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00457 00458 CALL sg_nuclear ( omat(1:m,1:m), l, basis%am(1:m,l), basis%am(1:m,l) ) 00459 CALL contract2(integrals%core(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00460 00461 DEALLOCATE (omat,STAT=ierr) 00462 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00463 END DO 00464 CASE (GTH_PSEUDO) 00465 alpha = 1._dp/potential%gth_pot%rc/SQRT(2._dp) 00466 DO l=0,3 00467 n = integrals%n(l) 00468 m = basis%nprim(l) 00469 IF(n>0 .AND. m>0) THEN 00470 ALLOCATE (omat(m,m),spmat(n,5),xmat(m),STAT=ierr) 00471 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00472 00473 omat = 0._dp 00474 CALL sg_erf ( omat(1:m,1:m), l, alpha, basis%am(1:m,l), basis%am(1:m,l) ) 00475 omat(1:m,1:m) = -potential%gth_pot%zion*omat(1:m,1:m) 00476 CALL contract2(integrals%core(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00477 DO i=1,potential%gth_pot%ncl 00478 omat = 0._dp 00479 CALL sg_gpot ( omat(1:m,1:m), i-1, potential%gth_pot%rc, l, basis%am(1:m,l), basis%am(1:m,l) ) 00480 omat(1:m,1:m) = potential%gth_pot%cl(i)*omat(1:m,1:m) 00481 CALL contract2add(integrals%core(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00482 END DO 00483 IF (potential%gth_pot%lpotextended) THEN 00484 DO k=1,potential%gth_pot%nexp_lpot 00485 DO i=1,potential%gth_pot%nct_lpot(k) 00486 omat = 0._dp 00487 CALL sg_gpot ( omat(1:m,1:m), i-1, potential%gth_pot%alpha_lpot(k), l, & 00488 basis%am(1:m,l), basis%am(1:m,l) ) 00489 omat(1:m,1:m) = potential%gth_pot%cval_lpot(i,k)*omat(1:m,1:m) 00490 CALL contract2add(integrals%core(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00491 END DO 00492 END DO 00493 END IF 00494 IF (potential%gth_pot%lsdpot) THEN 00495 DO k=1,potential%gth_pot%nexp_lsd 00496 DO i=1,potential%gth_pot%nct_lsd(k) 00497 omat = 0._dp 00498 CALL sg_gpot ( omat(1:m,1:m), i-1, potential%gth_pot%alpha_lsd(k), l, & 00499 basis%am(1:m,l), basis%am(1:m,l) ) 00500 omat(1:m,1:m) = potential%gth_pot%cval_lsd(i,k)*omat(1:m,1:m) 00501 CALL contract2add(integrals%clsd(1:n,1:n,l),omat(1:m,1:m),basis%cm(1:m,1:n,l), error) 00502 END DO 00503 END DO 00504 END IF 00505 00506 spmat = 0._dp 00507 k = potential%gth_pot%nl(l) 00508 DO i=1,k 00509 CALL sg_proj_ol ( xmat(1:m), l, basis%am(1:m,l), i-1, potential%gth_pot%rcnl(l) ) 00510 spmat(1:n,i) = MATMUL(TRANSPOSE(basis%cm(1:m,1:n,l)),xmat(1:m)) 00511 END DO 00512 IF(k>0) THEN 00513 integrals%hnl(1:n,1:n,l) = MATMUL(spmat(1:n,1:k),& 00514 MATMUL(potential%gth_pot%hnl(1:k,1:k,l),TRANSPOSE(spmat(1:n,1:k)))) 00515 END IF 00516 00517 DEALLOCATE (omat,spmat,xmat,STAT=ierr) 00518 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00519 END IF 00520 END DO 00521 CASE DEFAULT 00522 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00523 END SELECT 00524 00525 CASE (STO_BASIS) 00526 00527 SELECT CASE (potential%ppot_type) 00528 CASE (NO_PSEUDO) 00529 DO l=0,3 00530 n = integrals%n(l) 00531 CALL sto_nuclear ( integrals%core(1:n,1:n,l), basis%ns(1:n,l), basis%as(1:n,l),& 00532 basis%ns(1:n,l), basis%as(1:n,l) ) 00533 END DO 00534 CASE (GTH_PSEUDO) 00535 rad => basis%grid%rad 00536 m = basis%grid%nr 00537 ALLOCATE (cpot(1:m),STAT=ierr) 00538 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00539 rc = potential%gth_pot%rc 00540 alpha = 1._dp/rc/SQRT(2._dp) 00541 00542 ! local pseudopotential, we use erf = 1/r - erfc 00543 integrals%core = 0._dp 00544 DO i=1,m 00545 cpot(i) = potential%gth_pot%zion*erfc(alpha*rad(i))/rad(i) 00546 END DO 00547 DO i=1,potential%gth_pot%ncl 00548 ii = 2*(i-1) 00549 cpot(1:m) = cpot(1:m) + potential%gth_pot%cl(i) * (rad/rc)**ii * EXP(-0.5_dp*(rad/rc)**2) 00550 END DO 00551 IF (potential%gth_pot%lpotextended) THEN 00552 DO k=1,potential%gth_pot%nexp_lpot 00553 al = potential%gth_pot%alpha_lpot(k) 00554 DO i=1,potential%gth_pot%nct_lpot(k) 00555 ii = 2*(i-1) 00556 cpot(1:m) = cpot(1:m) + potential%gth_pot%cval_lpot(i,k)*(rad/al)**ii * EXP(-0.5_dp*(rad/al)**2) 00557 END DO 00558 END DO 00559 END IF 00560 CALL numpot_matrix(integrals%core,cpot,basis,0,error) 00561 DO l=0,3 00562 n = integrals%n(l) 00563 ALLOCATE (omat(n,n),STAT=ierr) 00564 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00565 omat = 0._dp 00566 CALL sto_nuclear ( omat(1:n,1:n), basis%ns(1:n,l), basis%as(1:n,l),& 00567 basis%ns(1:n,l), basis%as(1:n,l) ) 00568 integrals%core(1:n,1:n,l) = integrals%core(1:n,1:n,l) - potential%gth_pot%zion*omat(1:n,1:n) 00569 DEALLOCATE (omat,STAT=ierr) 00570 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00571 END DO 00572 00573 IF (potential%gth_pot%lsdpot) THEN 00574 cpot = 0._dp 00575 DO k=1,potential%gth_pot%nexp_lsd 00576 al = potential%gth_pot%alpha_lsd(k) 00577 DO i=1,potential%gth_pot%nct_lsd(k) 00578 ii = 2*(i-1) 00579 cpot = cpot + potential%gth_pot%cval_lsd(i,k)*(rad/al)**ii * EXP(-0.5_dp*(rad/al)**2) 00580 END DO 00581 END DO 00582 CALL numpot_matrix(integrals%clsd,cpot,basis,0,error) 00583 END IF 00584 00585 DO l=0,3 00586 n = integrals%n(l) 00587 ! non local pseudopotential 00588 ALLOCATE (spmat(n,5),STAT=ierr) 00589 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00590 spmat = 0._dp 00591 k = potential%gth_pot%nl(l) 00592 DO i=1,k 00593 rc = potential%gth_pot%rcnl(l) 00594 cpot = sqrt2/SQRT(gamma1(l+2*i-1))*rad**(l+2*i-2) * EXP(-0.5_dp*(rad/rc)**2) / rc**(l+2*i-0.5_dp) 00595 DO j=1,basis%nbas(l) 00596 spmat(j,i) = integrate_grid(cpot,basis%bf(:,j,l),basis%grid) 00597 END DO 00598 END DO 00599 integrals%hnl(1:n,1:n,l) = MATMUL(spmat(1:n,1:k),& 00600 MATMUL(potential%gth_pot%hnl(1:k,1:k,l),TRANSPOSE(spmat(1:n,1:k)))) 00601 00602 DEALLOCATE (spmat,STAT=ierr) 00603 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00604 00605 END DO 00606 00607 DEALLOCATE (cpot,STAT=ierr) 00608 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00609 00610 CASE DEFAULT 00611 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00612 END SELECT 00613 00614 CASE (NUM_BASIS) 00615 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00616 END SELECT 00617 00618 END IF 00619 00620 CALL timestop(handle) 00621 00622 END SUBROUTINE atom_ppint_setup 00623 ! ***************************************************************************** 00624 SUBROUTINE atom_relint_setup(integrals,basis,reltyp,zcore,error) 00625 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00626 TYPE(atom_basis_type), INTENT(INOUT) :: basis 00627 INTEGER, INTENT(IN) :: reltyp 00628 REAL(dp), OPTIONAL :: zcore 00629 TYPE(cp_error_type), INTENT(INOUT) :: error 00630 00631 CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_setup', 00632 routineP = moduleN//':'//routineN 00633 00634 INTEGER :: dkhorder, handle, ierr, l, m, 00635 n 00636 LOGICAL :: failure 00637 REAL(dp) :: cspeed 00638 REAL(dp), ALLOCATABLE, DIMENSION(:) :: cpot 00639 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: modpot 00640 REAL(KIND=dp), ALLOCATABLE, 00641 DIMENSION(:, :, :) :: pvp, sp, tp, vp 00642 00643 failure = .FALSE. 00644 00645 CALL timeset(routineN,handle) 00646 00647 cspeed = 1._dp/a_fine 00648 00649 SELECT CASE (reltyp) 00650 CASE DEFAULT 00651 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00652 CASE (do_nonrel_atom,do_zoramp_atom) 00653 dkhorder = -1 00654 CASE (do_dkh0_atom) 00655 dkhorder = 0 00656 CASE (do_dkh1_atom) 00657 dkhorder = 1 00658 CASE (do_dkh2_atom) 00659 dkhorder = 2 00660 CASE (do_dkh3_atom) 00661 dkhorder = 3 00662 CASE (do_dkh4_atom) 00663 dkhorder = 4 00664 CASE (do_dkh5_atom) 00665 dkhorder = 5 00666 END SELECT 00667 00668 SELECT CASE (reltyp) 00669 CASE DEFAULT 00670 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00671 CASE (do_nonrel_atom) 00672 ! nothing to do 00673 NULLIFY (integrals%tzora,integrals%hdkh) 00674 CASE (do_zoramp_atom) 00675 00676 NULLIFY (integrals%hdkh) 00677 00678 IF ( integrals%zorastat == 0 ) THEN 00679 n = MAXVAL(basis%nbas) 00680 ALLOCATE (integrals%tzora(n,n,0:3),STAT=ierr) 00681 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00682 integrals%tzora = 0._dp 00683 m = basis%grid%nr 00684 ALLOCATE (modpot(1:m),cpot(1:m),STAT=ierr) 00685 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00686 CALL calculate_model_potential(modpot,basis%grid,zcore,error) 00687 ! Zora potential 00688 cpot(1:m) = modpot(1:m)/(4._dp*cspeed*cspeed - 2._dp*modpot(1:m)) 00689 cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m) 00690 CALL numpot_matrix(integrals%tzora,cpot,basis,0,error) 00691 DO l=0,3 00692 integrals%tzora(:,:,l) = REAL(l*(l+1),dp) * integrals%tzora(:,:,l) 00693 END DO 00694 cpot(1:m) = cpot(1:m)*basis%grid%rad2(1:m) 00695 CALL numpot_matrix(integrals%tzora,cpot,basis,2,error) 00696 ! 00697 DEALLOCATE (modpot,cpot,STAT=ierr) 00698 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00699 00700 integrals%zorastat = 1 00701 00702 END IF 00703 00704 CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_dkh4_atom, do_dkh5_atom) 00705 00706 NULLIFY (integrals%tzora) 00707 00708 CPPostcondition(PRESENT(zcore), cp_failure_level, routineP, error, failure) 00709 IF ( integrals%dkhstat == 0 ) THEN 00710 n = MAXVAL(basis%nbas) 00711 ALLOCATE (integrals%hdkh(n,n,0:3),STAT=ierr) 00712 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00713 integrals%hdkh = 0._dp 00714 00715 m = MAXVAL(basis%nprim) 00716 ALLOCATE (tp(m,m,0:3),sp(m,m,0:3),vp(m,m,0:3),pvp(m,m,0:3),STAT=ierr) 00717 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00718 tp = 0._dp 00719 sp = 0._dp 00720 vp = 0._dp 00721 pvp = 0._dp 00722 00723 SELECT CASE (basis%basis_type) 00724 CASE DEFAULT 00725 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00726 CASE (GTO_BASIS, CGTO_BASIS) 00727 00728 DO l=0,3 00729 m = basis%nprim(l) 00730 IF ( m > 0 ) THEN 00731 CALL sg_kinetic ( tp(1:m,1:m,l), l, basis%am(1:m,l), basis%am(1:m,l) ) 00732 CALL sg_overlap ( sp(1:m,1:m,l), l, basis%am(1:m,l), basis%am(1:m,l) ) 00733 CALL sg_nuclear ( vp(1:m,1:m,l), l, basis%am(1:m,l), basis%am(1:m,l) ) 00734 CALL sg_kinnuc ( pvp(1:m,1:m,l), l, basis%am(1:m,l), basis%am(1:m,l) ) 00735 vp(1:m,1:m,l) = -zcore*vp(1:m,1:m,l) 00736 pvp(1:m,1:m,l) = -zcore*pvp(1:m,1:m,l) 00737 END IF 00738 END DO 00739 00740 CASE (STO_BASIS) 00741 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00742 CASE (NUM_BASIS) 00743 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00744 END SELECT 00745 00746 CALL dkh_integrals(integrals,basis,dkhorder,sp,tp,vp,pvp,error) 00747 00748 integrals%dkhstat = 1 00749 00750 DEALLOCATE (tp,sp,vp,pvp,STAT=ierr) 00751 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00752 00753 ELSE 00754 CPPostcondition(ASSOCIATED(integrals%hdkh), cp_failure_level, routineP, error, failure) 00755 END IF 00756 00757 END SELECT 00758 00759 CALL timestop(handle) 00760 00761 END SUBROUTINE atom_relint_setup 00762 ! ***************************************************************************** 00763 SUBROUTINE dkh_integrals(integrals,basis,order,sp,tp,vp,pvp,error) 00764 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00765 TYPE(atom_basis_type), INTENT(INOUT) :: basis 00766 INTEGER, INTENT(IN) :: order 00767 REAL(dp), DIMENSION(:, :, 0:) :: sp, tp, vp, pvp 00768 TYPE(cp_error_type), INTENT(INOUT) :: error 00769 00770 CHARACTER(len=*), PARAMETER :: routineN = 'dkh_integrals', 00771 routineP = moduleN//':'//routineN 00772 00773 INTEGER :: l, m, n 00774 LOGICAL :: failure 00775 REAL(dp), DIMENSION(:, :, :), POINTER :: hdkh 00776 00777 failure = .FALSE. 00778 CPPrecondition(order>=0, cp_failure_level, routineP, error, failure) 00779 00780 hdkh => integrals%hdkh 00781 00782 DO l=0,3 00783 n = integrals%n(l) 00784 m = basis%nprim(l) 00785 IF ( n > 0 ) THEN 00786 CALL dkh_atom_transformation(sp(1:m,1:m,l),vp(1:m,1:m,l),tp(1:m,1:m,l),pvp(1:m,1:m,l),m,order) 00787 SELECT CASE (basis%basis_type) 00788 CASE DEFAULT 00789 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00790 CASE (GTO_BASIS) 00791 CPAssert(n==m,cp_failure_level,routineP,error,failure) 00792 integrals%hdkh(1:n,1:n,l) = tp(1:n,1:n,l) + vp(1:n,1:n,l) 00793 CASE (CGTO_BASIS) 00794 CALL contract2(integrals%hdkh(1:n,1:n,l),tp(1:m,1:m,l),basis%cm(1:m,1:n,l), error) 00795 CALL contract2add(integrals%hdkh(1:n,1:n,l),vp(1:m,1:m,l),basis%cm(1:m,1:n,l), error) 00796 CASE (STO_BASIS) 00797 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00798 CASE (NUM_BASIS) 00799 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00800 END SELECT 00801 ELSE 00802 integrals%hdkh(1:n,1:n,l) = 0._dp 00803 END IF 00804 END DO 00805 00806 END SUBROUTINE dkh_integrals 00807 ! ***************************************************************************** 00808 SUBROUTINE atom_int_release(integrals,error) 00809 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00810 TYPE(cp_error_type), INTENT(INOUT) :: error 00811 00812 CHARACTER(len=*), PARAMETER :: routineN = 'atom_int_release', 00813 routineP = moduleN//':'//routineN 00814 00815 INTEGER :: ierr, ll 00816 LOGICAL :: failure 00817 00818 failure = .FALSE. 00819 00820 IF ( ASSOCIATED(integrals%ovlp) ) THEN 00821 DEALLOCATE (integrals%ovlp,STAT=ierr) 00822 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00823 END IF 00824 IF ( ASSOCIATED(integrals%kin) ) THEN 00825 DEALLOCATE (integrals%kin,STAT=ierr) 00826 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00827 END IF 00828 IF ( ASSOCIATED(integrals%conf) ) THEN 00829 DEALLOCATE (integrals%conf,STAT=ierr) 00830 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00831 END IF 00832 DO ll=1,SIZE(integrals%ceri) 00833 IF ( ASSOCIATED(integrals%ceri(ll)%int) ) THEN 00834 DEALLOCATE (integrals%ceri(ll)%int,STAT=ierr) 00835 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00836 END IF 00837 IF ( ASSOCIATED(integrals%eeri(ll)%int) ) THEN 00838 DEALLOCATE (integrals%eeri(ll)%int,STAT=ierr) 00839 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00840 END IF 00841 END DO 00842 IF ( ASSOCIATED(integrals%utrans) ) THEN 00843 DEALLOCATE (integrals%utrans,STAT=ierr) 00844 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00845 END IF 00846 IF ( ASSOCIATED(integrals%uptrans) ) THEN 00847 DEALLOCATE (integrals%uptrans,STAT=ierr) 00848 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00849 END IF 00850 00851 integrals%status = 0 00852 00853 END SUBROUTINE atom_int_release 00854 ! ***************************************************************************** 00855 SUBROUTINE atom_ppint_release(integrals,error) 00856 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00857 TYPE(cp_error_type), INTENT(INOUT) :: error 00858 00859 CHARACTER(len=*), PARAMETER :: routineN = 'atom_ppint_release', 00860 routineP = moduleN//':'//routineN 00861 00862 INTEGER :: ierr 00863 LOGICAL :: failure 00864 00865 failure = .FALSE. 00866 00867 IF ( ASSOCIATED(integrals%hnl) ) THEN 00868 DEALLOCATE (integrals%hnl,STAT=ierr) 00869 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00870 END IF 00871 IF ( ASSOCIATED(integrals%core) ) THEN 00872 DEALLOCATE (integrals%core,STAT=ierr) 00873 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00874 END IF 00875 IF ( ASSOCIATED(integrals%clsd) ) THEN 00876 DEALLOCATE (integrals%clsd,STAT=ierr) 00877 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00878 END IF 00879 00880 integrals%ppstat = 0 00881 00882 END SUBROUTINE atom_ppint_release 00883 ! ***************************************************************************** 00884 SUBROUTINE atom_relint_release(integrals,error) 00885 TYPE(atom_integrals), INTENT(INOUT) :: integrals 00886 TYPE(cp_error_type), INTENT(INOUT) :: error 00887 00888 CHARACTER(len=*), PARAMETER :: routineN = 'atom_relint_release', 00889 routineP = moduleN//':'//routineN 00890 00891 INTEGER :: ierr 00892 LOGICAL :: failure 00893 00894 failure = .FALSE. 00895 00896 IF ( ASSOCIATED(integrals%tzora) ) THEN 00897 DEALLOCATE (integrals%tzora,STAT=ierr) 00898 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00899 END IF 00900 IF ( ASSOCIATED(integrals%hdkh) ) THEN 00901 DEALLOCATE (integrals%hdkh,STAT=ierr) 00902 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00903 END IF 00904 00905 integrals%zorastat = 0 00906 integrals%dkhstat = 0 00907 00908 END SUBROUTINE atom_relint_release 00909 ! ***************************************************************************** 00910 SUBROUTINE calculate_model_potential(modpot,grid,zcore,error) 00911 REAL(dp), DIMENSION(:), INTENT(INOUT) :: modpot 00912 TYPE(grid_atom_type), INTENT(IN) :: grid 00913 REAL(dp), INTENT(IN) :: zcore 00914 TYPE(cp_error_type), INTENT(INOUT) :: error 00915 00916 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_model_potential', 00917 routineP = moduleN//':'//routineN 00918 00919 INTEGER :: i, ierr, ii, l, ll, n, z 00920 LOGICAL :: failure 00921 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pot, rho 00922 TYPE(atom_state) :: state 00923 00924 n = SIZE(modpot) 00925 ALLOCATE(rho(n),pot(n),STAT=ierr) 00926 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00927 00928 ! fill default occupation 00929 state%core = 0._dp 00930 state%occ = 0._dp 00931 state%multiplicity = -1 00932 z = NINT(zcore) 00933 DO l=0,3 00934 IF ( ptable(z)%e_conv(l) /= 0 ) THEN 00935 state%maxl_occ = l 00936 ll = 2*(2*l+1) 00937 DO i=1,SIZE(state%occ,2) 00938 ii = ptable(z)%e_conv(l) - (i-1)*ll 00939 IF ( ii <= ll ) THEN 00940 state%occ(l,i) = ii 00941 EXIT 00942 ELSE 00943 state%occ(l,i) = ll 00944 END IF 00945 END DO 00946 END IF 00947 END DO 00948 00949 modpot = -zcore/grid%rad(:) 00950 00951 ! Coulomb potential 00952 CALL slater_density(rho,pot,NINT(zcore),state,grid,error) 00953 CALL coulomb_potential_numeric(pot,rho,grid,error) 00954 modpot = modpot + pot 00955 00956 ! XC potential 00957 CALL wigner_slater_functional(rho,pot,error) 00958 modpot = modpot + pot 00959 00960 DEALLOCATE(rho,pot,STAT=ierr) 00961 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00962 00963 END SUBROUTINE calculate_model_potential 00964 ! ***************************************************************************** 00965 SUBROUTINE contract2 ( int, omat, cm, error ) 00966 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int 00967 REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm 00968 TYPE(cp_error_type), INTENT(INOUT) :: error 00969 00970 CHARACTER(len=*), PARAMETER :: routineN = 'contract2', 00971 routineP = moduleN//':'//routineN 00972 00973 INTEGER :: handle, m, n 00974 00975 CALL timeset(routineN,handle) 00976 00977 n = SIZE(int,1) 00978 m = SIZE(omat,1) 00979 00980 INT(1:n,1:n) = MATMUL(TRANSPOSE(cm(1:m,1:n)),MATMUL(omat(1:m,1:m),cm(1:m,1:n))) 00981 00982 CALL timestop(handle) 00983 00984 END SUBROUTINE contract2 00985 ! ***************************************************************************** 00986 SUBROUTINE contract2add ( int, omat, cm, error ) 00987 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int 00988 REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm 00989 TYPE(cp_error_type), INTENT(INOUT) :: error 00990 00991 CHARACTER(len=*), PARAMETER :: routineN = 'contract2add', 00992 routineP = moduleN//':'//routineN 00993 00994 INTEGER :: handle, m, n 00995 00996 CALL timeset(routineN,handle) 00997 00998 n = SIZE(int,1) 00999 m = SIZE(omat,1) 01000 01001 INT(1:n,1:n) = INT(1:n,1:n) + MATMUL(TRANSPOSE(cm(1:m,1:n)),MATMUL(omat(1:m,1:m),cm(1:m,1:n))) 01002 01003 CALL timestop(handle) 01004 01005 END SUBROUTINE contract2add 01006 ! ***************************************************************************** 01007 SUBROUTINE contract4 ( eri, omat, cm1, cm2, error ) 01008 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: eri 01009 REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm1, cm2 01010 TYPE(cp_error_type), INTENT(INOUT) :: error 01011 01012 CHARACTER(len=*), PARAMETER :: routineN = 'contract4', 01013 routineP = moduleN//':'//routineN 01014 01015 INTEGER :: handle, i1, i2, ierr, m1, m2, 01016 mm1, mm2, n1, n2, nn1, nn2 01017 LOGICAL :: failure = .FALSE. 01018 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: amat, atran, bmat, btran, hint 01019 01020 CALL timeset(routineN,handle) 01021 01022 m1 = SIZE(cm1,1) 01023 n1 = SIZE(cm1,2) 01024 m2 = SIZE(cm2,1) 01025 n2 = SIZE(cm2,2) 01026 nn1 = SIZE(eri,1) 01027 nn2 = SIZE(eri,2) 01028 mm1 = SIZE(omat,1) 01029 mm2 = SIZE(omat,2) 01030 01031 ALLOCATE(amat(m1,m1),atran(n1,n1),bmat(m2,m2),btran(n2,n2),STAT=ierr) 01032 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 01033 ALLOCATE(hint(mm1,nn2),STAT=ierr) 01034 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 01035 01036 DO i1=1,mm1 01037 CALL iunpack(bmat(1:m2,1:m2),omat(i1,1:mm2),m2) 01038 CALL contract2( btran(1:n2,1:n2), bmat(1:m2,1:m2), cm2(1:m2,1:n2), error ) 01039 CALL ipack(btran(1:n2,1:n2),hint(i1,1:nn2),n2) 01040 END DO 01041 01042 DO i2=1,nn2 01043 CALL iunpack(amat(1:m1,1:m1),hint(1:mm1,i2),m1) 01044 CALL contract2( atran(1:n1,1:n1), amat(1:m1,1:m1), cm1(1:m1,1:n1), error ) 01045 CALL ipack(atran(1:n1,1:n1),eri(1:nn1,i2),n1) 01046 END DO 01047 01048 DEALLOCATE(amat,atran,bmat,btran,STAT=ierr) 01049 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 01050 DEALLOCATE(hint,STAT=ierr) 01051 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 01052 01053 CALL timestop(handle) 01054 01055 END SUBROUTINE contract4 01056 ! ***************************************************************************** 01057 SUBROUTINE ipack(mat,vec,n) 01058 REAL(dp), DIMENSION(:, :), INTENT(IN) :: mat 01059 REAL(dp), DIMENSION(:), INTENT(INOUT) :: vec 01060 INTEGER, INTENT(IN) :: n 01061 01062 INTEGER :: i, ij, j 01063 01064 ij = 0 01065 DO i=1,n 01066 DO j=i,n 01067 ij = ij + 1 01068 vec(ij) = mat(i,j) 01069 END DO 01070 END DO 01071 01072 END SUBROUTINE ipack 01073 01074 SUBROUTINE iunpack(mat,vec,n) 01075 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: mat 01076 REAL(dp), DIMENSION(:), INTENT(IN) :: vec 01077 INTEGER, INTENT(IN) :: n 01078 01079 INTEGER :: i, ij, j 01080 01081 ij = 0 01082 DO i=1,n 01083 DO j=i,n 01084 ij = ij + 1 01085 mat(i,j) = vec(ij) 01086 mat(j,i) = vec(ij) 01087 END DO 01088 END DO 01089 01090 END SUBROUTINE iunpack 01091 ! ***************************************************************************** 01092 01093 END MODULE atom_operators
1.7.3