|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00010 MODULE qs_dftb_matrices 00011 USE array_types, ONLY: array_i1d_obj,& 00012 array_new,& 00013 array_nullify,& 00014 array_release 00015 USE atomic_kind_types, ONLY: atomic_kind_type,& 00016 get_atomic_kind,& 00017 get_atomic_kind_set,& 00018 is_hydrogen 00019 USE atprop_types, ONLY: atprop_array_init 00020 USE block_p_types, ONLY: block_p_type 00021 USE cp_control_types, ONLY: dft_control_type,& 00022 dftb_control_type 00023 USE cp_dbcsr_interface, ONLY: cp_dbcsr_add,& 00024 cp_dbcsr_copy,& 00025 cp_dbcsr_create,& 00026 cp_dbcsr_finalize,& 00027 cp_dbcsr_get_block_p,& 00028 cp_dbcsr_init,& 00029 cp_dbcsr_multiply 00030 USE cp_dbcsr_operations, ONLY: cp_dbcsr_add_block_node,& 00031 cp_dbcsr_allocate_matrix_set 00032 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 00033 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00034 cp_dbcsr_type 00035 USE cp_output_handling, ONLY: cp_p_file,& 00036 cp_print_key_finished_output,& 00037 cp_print_key_should_output,& 00038 cp_print_key_unit_nr 00039 USE cp_para_types, ONLY: cp_para_env_type 00040 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00041 dbcsr_type_antisymmetric,& 00042 dbcsr_type_symmetric 00043 USE dbcsr_util, ONLY: convert_offsets_to_sizes 00044 USE f77_blas 00045 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00046 section_vals_type 00047 USE kinds, ONLY: dp 00048 USE message_passing, ONLY: mp_sum 00049 USE mulliken, ONLY: mulliken_charges 00050 USE particle_types, ONLY: get_particle_set,& 00051 particle_type 00052 USE qs_dftb_coulomb, ONLY: build_dftb_coulomb 00053 USE qs_dftb_types, ONLY: qs_dftb_atom_type,& 00054 qs_dftb_pairpot_type 00055 USE qs_dftb_utils, ONLY: get_dftb_atom_param 00056 USE qs_energy_types, ONLY: qs_energy_type 00057 USE qs_environment_types, ONLY: get_qs_env,& 00058 qs_environment_type,& 00059 set_qs_env 00060 USE qs_force_types, ONLY: qs_force_type 00061 USE qs_ks_types, ONLY: qs_ks_env_type 00062 USE qs_mo_types, ONLY: get_mo_set,& 00063 mo_set_p_type 00064 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 00065 neighbor_list_iterate,& 00066 neighbor_list_iterator_create,& 00067 neighbor_list_iterator_p_type,& 00068 neighbor_list_iterator_release,& 00069 neighbor_list_set_p_type 00070 USE qs_rho_types, ONLY: qs_rho_type 00071 USE timings, ONLY: timeset,& 00072 timestop 00073 USE virial_methods, ONLY: virial_pair_force 00074 USE virial_types, ONLY: virial_type 00075 #include "cp_common_uses.h" 00076 00077 IMPLICIT NONE 00078 00079 INTEGER,DIMENSION(16),PARAMETER :: orbptr = (/ 0, 1, 1, 1, 00080 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 /) 00081 00082 ! Maximum number of points used for interpolation 00083 INTEGER, PARAMETER :: max_inter = 5 00084 ! Maximum number of points used for extrapolation 00085 INTEGER, PARAMETER :: max_extra = 9 00086 ! see also qs_dftb_parameters 00087 REAL(dp), PARAMETER :: slako_d0 = 1._dp 00088 ! pointer to skab 00089 INTEGER, DIMENSION(0:3,0:3,0:3,0:3,0:3):: iptr 00090 ! screening for gamma function 00091 REAL(dp), PARAMETER :: tol_gamma = 1.e-4_dp 00092 ! small real number 00093 REAL(dp), PARAMETER :: rtiny = 1.e-10_dp 00094 00095 PRIVATE 00096 00097 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_matrices' 00098 00099 PUBLIC :: build_dftb_matrices, build_dftb_ks_matrix 00100 00101 CONTAINS 00102 00103 ! ***************************************************************************** 00104 SUBROUTINE build_dftb_matrices(qs_env,para_env,calculate_forces,error) 00105 00106 TYPE(qs_environment_type), POINTER :: qs_env 00107 TYPE(cp_para_env_type), POINTER :: para_env 00108 LOGICAL, INTENT(IN) :: calculate_forces 00109 TYPE(cp_error_type), INTENT(inout) :: error 00110 00111 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dftb_matrices', 00112 routineP = moduleN//':'//routineN 00113 00114 INTEGER :: atom_a, atom_b, handle, i, iatom, icol, ikind, inode, irow, 00115 istat, iw, jatom, jkind, l1, l2, la, last_jatom, lb, llm, lmaxi, lmaxj, 00116 m, n1, n2, n_urpoly, natom, natorb_a, natorb_b, ncol, neighbor_list_id, 00117 ngrd, ngrdcut, nkind, nmat, nrow, spdim 00118 INTEGER, DIMENSION(:), POINTER :: atom_of_kind, felem, lelem, 00119 rbs 00120 LOGICAL :: defined, failure, found, 00121 hb_sr_damp, new_atom_b, 00122 use_virial 00123 REAL(KIND=dp) :: ddr, dgam, dgrd, dr, drm, 00124 drp, erep, erepij, f0, f1, 00125 foab, fow, ga, gb, s_cut, 00126 urep_cut 00127 REAL(KIND=dp), DIMENSION(0:3) :: eta_a, eta_b, skself 00128 REAL(KIND=dp), DIMENSION(10) :: urep 00129 REAL(KIND=dp), DIMENSION(2) :: surr 00130 REAL(KIND=dp), DIMENSION(20) :: skabij, skabji 00131 REAL(KIND=dp), DIMENSION(3) :: force_ab, force_rr, force_w, 00132 rij, srep 00133 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dfblock, dsblock, fblock, 00134 fmatij, fmatji, gblock, pblock, sblock, scoeff, smatij, smatji, spxr, 00135 wblock 00136 TYPE(array_i1d_obj) :: row_blk_sizes 00137 TYPE(atomic_kind_type), DIMENSION(:), 00138 POINTER :: atomic_kind_set 00139 TYPE(atomic_kind_type), POINTER :: atomic_kind 00140 TYPE(block_p_type), DIMENSION(2:4) :: dgblocks, dsblocks 00141 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00142 POINTER :: gamma_matrix, matrix_h, 00143 matrix_p, matrix_s, matrix_w 00144 TYPE(cp_logger_type), POINTER :: logger 00145 TYPE(dbcsr_distribution_obj), POINTER :: dbcsr_dist 00146 TYPE(dft_control_type), POINTER :: dft_control 00147 TYPE(dftb_control_type), POINTER :: dftb_control 00148 TYPE(neighbor_list_iterator_p_type), 00149 DIMENSION(:), POINTER :: nl_iterator 00150 TYPE(neighbor_list_set_p_type), 00151 DIMENSION(:), POINTER :: sab_orb 00152 TYPE(particle_type), DIMENSION(:), 00153 POINTER :: particle_set 00154 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b 00155 TYPE(qs_dftb_pairpot_type), 00156 DIMENSION(:, :), POINTER :: dftb_potential 00157 TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji 00158 TYPE(qs_energy_type), POINTER :: energy 00159 TYPE(qs_force_type), DIMENSION(:), 00160 POINTER :: force 00161 TYPE(qs_rho_type), POINTER :: rho 00162 TYPE(virial_type), POINTER :: virial 00163 00164 CALL timeset(routineN,handle) 00165 00166 ! set pointers 00167 iptr = 0 00168 DO la=0,3 00169 DO lb=0,3 00170 llm=0 00171 DO l1=0,MAX(la,lb) 00172 DO l2=0,MIN(l1,la,lb) 00173 DO m=0,l2 00174 llm=llm+1 00175 iptr(l1,l2,m,la,lb)=llm 00176 END DO 00177 END DO 00178 END DO 00179 END DO 00180 END DO 00181 00182 NULLIFY(logger) 00183 logger => cp_error_get_logger(error) 00184 00185 ! Allocate the overlap and Hamiltonian matrix 00186 CALL setup_matrices(qs_env,calculate_forces,error) 00187 00188 NULLIFY (matrix_h,matrix_s,matrix_p,matrix_w,gamma_matrix,atomic_kind_set,sab_orb) 00189 00190 CALL get_qs_env(qs_env=qs_env,& 00191 energy=energy,& 00192 atomic_kind_set=atomic_kind_set,& 00193 matrix_h=matrix_h,matrix_s=matrix_s,& 00194 dft_control=dft_control,error=error) 00195 00196 dftb_control => dft_control%qs_control%dftb_control 00197 00198 NULLIFY (dftb_potential) 00199 CALL get_qs_env(qs_env=qs_env,& 00200 dftb_potential=dftb_potential,error=error) 00201 00202 ! gamma matrix allocation 00203 IF ( dftb_control%self_consistent ) THEN 00204 IF(calculate_forces) THEN 00205 nmat=4 00206 ELSE 00207 nmat=1 00208 END IF 00209 CALL get_qs_env(qs_env=qs_env,& 00210 particle_set=particle_set,& 00211 neighbor_list_id=neighbor_list_id,& 00212 dbcsr_dist=dbcsr_dist,error=error) 00213 natom = SIZE(particle_set) 00214 nrow = natom 00215 ncol = natom 00216 ALLOCATE (felem(natom),lelem(natom),STAT=istat) 00217 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00218 DO iatom = 1, natom 00219 felem(iatom) = iatom 00220 lelem(iatom) = iatom 00221 ENDDO 00222 CALL get_qs_env(qs_env=qs_env,& 00223 gamma_matrix=gamma_matrix,error=error) 00224 00225 ALLOCATE (rbs(natom), STAT=istat) 00226 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00227 CALL convert_offsets_to_sizes (felem, rbs, lelem) 00228 CALL array_nullify (row_blk_sizes) 00229 CALL array_new (row_blk_sizes, rbs, gift=.TRUE.) 00230 00231 00232 CALL cp_dbcsr_allocate_matrix_set(gamma_matrix,nmat,error=error) 00233 ALLOCATE(gamma_matrix(1)%matrix) 00234 CALL cp_dbcsr_init(gamma_matrix(1)%matrix, error=error) 00235 00236 CALL cp_dbcsr_create(matrix=gamma_matrix(1)%matrix, & 00237 name="GAMMA MATRIX", & 00238 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric,& 00239 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00240 nblks=0, nze=0, mutable_work=.TRUE., & 00241 error=error) 00242 00243 DO i=2,nmat 00244 ALLOCATE(gamma_matrix(i)%matrix) 00245 CALL cp_dbcsr_init(gamma_matrix(i)%matrix, error=error) 00246 00247 CALL cp_dbcsr_create(matrix=gamma_matrix(i)%matrix, & 00248 name="DERIVATIVE GAMMA MATRIX", & 00249 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric,& 00250 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00251 nblks=0, nze=0, mutable_work=.TRUE., & 00252 error=error) 00253 END DO 00254 00255 CALL array_release (row_blk_sizes) 00256 00257 DEALLOCATE (felem,lelem,STAT=istat) 00258 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00259 END IF 00260 00261 IF(calculate_forces) THEN 00262 NULLIFY (rho,force,particle_set,matrix_w) 00263 CALL get_qs_env(qs_env=qs_env,& 00264 particle_set=particle_set,& 00265 rho=rho,matrix_w=matrix_w,& 00266 virial=virial,& 00267 force=force,error=error) 00268 matrix_p => rho%rho_ao 00269 00270 IF (SIZE(matrix_p) == 2) THEN 00271 CALL cp_dbcsr_add(matrix_p(1)%matrix,matrix_p(2)%matrix,& 00272 alpha_scalar=1.0_dp,beta_scalar=1.0_dp,error=error) 00273 CALL cp_dbcsr_add(matrix_w(1)%matrix,matrix_w(2)%matrix,& 00274 alpha_scalar=1.0_dp,beta_scalar=1.0_dp,error=error) 00275 END IF 00276 natom = SIZE (particle_set) 00277 ALLOCATE (atom_of_kind(natom),STAT=istat) 00278 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00279 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00280 atom_of_kind=atom_of_kind) 00281 use_virial = virial%pv_availability.AND.(.NOT.virial%pv_numer) 00282 END IF 00283 ! atomic energy decomposition 00284 IF (qs_env%atprop%energy) THEN 00285 CALL get_qs_env(qs_env=qs_env,particle_set=particle_set,error=error) 00286 natom = SIZE (particle_set) 00287 CALL atprop_array_init(qs_env%atprop%atecc,natom,error) 00288 END IF 00289 00290 erep = 0._dp 00291 00292 CALL get_qs_env(qs_env=qs_env,sab_orb=sab_orb,error=error) 00293 00294 nkind = SIZE(atomic_kind_set) 00295 00296 CALL neighbor_list_iterator_create(nl_iterator,sab_orb) 00297 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 00298 CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,& 00299 iatom=iatom,jatom=jatom,inode=inode,r=rij) 00300 00301 atomic_kind => atomic_kind_set(ikind) 00302 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00303 natom=natom,& 00304 dftb_parameter=dftb_kind_a) 00305 CALL get_dftb_atom_param(dftb_kind_a,& 00306 defined=defined,lmax=lmaxi,skself=skself,& 00307 eta=eta_a,natorb=natorb_a) 00308 00309 IF (.NOT.defined .OR. natorb_a < 1) CYCLE 00310 00311 atomic_kind => atomic_kind_set(jkind) 00312 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00313 dftb_parameter=dftb_kind_b) 00314 CALL get_dftb_atom_param(dftb_kind_b,& 00315 defined=defined,lmax=lmaxj,eta=eta_b,natorb=natorb_b) 00316 00317 IF (.NOT.defined .OR. natorb_b < 1) CYCLE 00318 00319 ! retrieve information on F and S matrix 00320 dftb_param_ij => dftb_potential(ikind,jkind) 00321 dftb_param_ji => dftb_potential(jkind,ikind) 00322 ! assume table size and type is symmetric 00323 ngrd = dftb_param_ij%ngrd 00324 ngrdcut = dftb_param_ij%ngrdcut 00325 dgrd = dftb_param_ij%dgrd 00326 ddr = dgrd*0.1_dp 00327 CPPrecondition(dftb_param_ij%llm==dftb_param_ji%llm,cp_failure_level,routineP,error,failure) 00328 llm = dftb_param_ij%llm 00329 fmatij => dftb_param_ij%fmat 00330 smatij => dftb_param_ij%smat 00331 fmatji => dftb_param_ji%fmat 00332 smatji => dftb_param_ji%smat 00333 ! repulsive pair potential 00334 n_urpoly = dftb_param_ij%n_urpoly 00335 urep_cut = dftb_param_ij%urep_cut 00336 urep = dftb_param_ij%urep 00337 spxr => dftb_param_ij%spxr 00338 scoeff => dftb_param_ij%scoeff 00339 spdim = dftb_param_ij%spdim 00340 s_cut = dftb_param_ij%s_cut 00341 srep = dftb_param_ij%srep 00342 surr = dftb_param_ij%surr 00343 00344 IF (inode==1) last_jatom=0 00345 00346 dr = SQRT(SUM(rij(:)**2)) 00347 IF (NINT(dr/dgrd) <= ngrdcut) THEN 00348 00349 IF (jatom /= last_jatom) THEN 00350 new_atom_b = .TRUE. 00351 last_jatom = jatom 00352 ELSE 00353 new_atom_b = .FALSE. 00354 END IF 00355 IF (new_atom_b) THEN 00356 icol = MAX(iatom,jatom) 00357 irow = MIN(iatom,jatom) 00358 NULLIFY (sblock) 00359 CALL cp_dbcsr_add_block_node(matrix=matrix_s(1)%matrix,& 00360 block_row=irow,& 00361 block_col=icol,& 00362 BLOCK=sblock,error=error) 00363 NULLIFY (fblock) 00364 CALL cp_dbcsr_add_block_node(matrix=matrix_h(1)%matrix,& 00365 block_row=irow,& 00366 block_col=icol,& 00367 BLOCK=fblock,error=error) 00368 IF ( dftb_control%self_consistent ) THEN 00369 NULLIFY (gblock) 00370 CALL cp_dbcsr_add_block_node(matrix=gamma_matrix(1)%matrix,& 00371 block_row=irow,& 00372 block_col=icol,& 00373 BLOCK=gblock,error=error) 00374 END IF 00375 IF (calculate_forces) THEN 00376 NULLIFY (pblock) 00377 CALL cp_dbcsr_get_block_p(matrix=matrix_p(1)%matrix,& 00378 row=irow,col=icol,block=pblock,found=found) 00379 CPPrecondition(ASSOCIATED(pblock),cp_failure_level,routineP,error,failure) 00380 NULLIFY (wblock) 00381 CALL cp_dbcsr_get_block_p(matrix=matrix_w(1)%matrix,& 00382 row=irow,col=icol,block=wblock,found=found) 00383 CPPrecondition(ASSOCIATED(wblock),cp_failure_level,routineP,error,failure) 00384 IF ( dftb_control%self_consistent ) THEN 00385 DO i=2,4 00386 NULLIFY(dsblocks(i)%block) 00387 CALL cp_dbcsr_add_block_node(matrix=matrix_s(i)%matrix,& 00388 block_row=irow,& 00389 block_col=icol,& 00390 BLOCK=dsblocks(i)%block,error=error) 00391 NULLIFY (dgblocks(i)%block) 00392 CALL cp_dbcsr_add_block_node(matrix=gamma_matrix(i)%matrix,& 00393 block_row=irow,& 00394 block_col=icol,& 00395 BLOCK=dgblocks(i)%block,error=error) 00396 END DO 00397 END IF 00398 END IF 00399 END IF 00400 00401 IF (iatom == jatom .AND. dr < 0.001_dp) THEN 00402 ! diagonal block 00403 DO i=1,natorb_a 00404 sblock(i,i) = sblock(i,i) + 1._dp 00405 fblock(i,i) = fblock(i,i) + skself(orbptr(i)) 00406 END DO 00407 ELSE 00408 ! off-diagonal block 00409 IF ( irow == iatom ) THEN 00410 CALL getskz(smatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00411 CALL getskz(smatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00412 CALL turnsk(sblock,skabji,skabij,rij,dr,lmaxi,lmaxj) 00413 CALL getskz(fmatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00414 CALL getskz(fmatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00415 CALL turnsk(fblock,skabji,skabij,rij,dr,lmaxi,lmaxj) 00416 ELSE 00417 CALL getskz(smatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00418 CALL getskz(smatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00419 CALL turnsk(sblock,skabij,skabji,-rij,dr,lmaxj,lmaxi) 00420 CALL getskz(fmatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00421 CALL getskz(fmatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00422 CALL turnsk(fblock,skabij,skabji,-rij,dr,lmaxj,lmaxi) 00423 END IF 00424 IF(calculate_forces) THEN 00425 force_ab = 0._dp 00426 force_w = 0._dp 00427 n1 = SIZE(fblock,1) 00428 n2 = SIZE(fblock,2) 00429 ALLOCATE (dfblock(n1,n2),dsblock(n1,n2),STAT=istat) 00430 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00431 DO i=1,3 00432 rij(i) = rij(i) + ddr 00433 dr = SQRT(SUM(rij(:)**2)) 00434 dfblock=0._dp 00435 dsblock=0._dp 00436 IF ( irow == iatom ) THEN 00437 CALL getskz(smatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00438 CALL getskz(smatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00439 CALL turnsk(dsblock,skabji,skabij,rij,dr,lmaxi,lmaxj) 00440 CALL getskz(fmatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00441 CALL getskz(fmatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00442 CALL turnsk(dfblock,skabji,skabij,rij,dr,lmaxi,lmaxj) 00443 ELSE 00444 CALL getskz(smatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00445 CALL getskz(smatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00446 CALL turnsk(dsblock,skabij,skabji,-rij,dr,lmaxj,lmaxi) 00447 CALL getskz(fmatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00448 CALL getskz(fmatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00449 CALL turnsk(dfblock,skabij,skabji,-rij,dr,lmaxj,lmaxi) 00450 END IF 00451 rij(i) = rij(i) - 2._dp*ddr 00452 dr = SQRT(SUM(rij(:)**2)) 00453 dsblock = -dsblock 00454 dfblock = -dfblock 00455 IF ( irow == iatom ) THEN 00456 CALL getskz(smatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00457 CALL getskz(smatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00458 CALL turnsk(dsblock,skabji,skabij,rij,dr,lmaxi,lmaxj) 00459 CALL getskz(fmatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00460 CALL getskz(fmatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00461 CALL turnsk(dfblock,skabji,skabij,rij,dr,lmaxi,lmaxj) 00462 ELSE 00463 CALL getskz(smatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00464 CALL getskz(smatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00465 CALL turnsk(dsblock,skabij,skabji,-rij,dr,lmaxj,lmaxi) 00466 CALL getskz(fmatij,skabij,dr,ngrd,ngrdcut,dgrd,llm) 00467 CALL getskz(fmatji,skabji,dr,ngrd,ngrdcut,dgrd,llm) 00468 CALL turnsk(dfblock,skabij,skabji,-rij,dr,lmaxj,lmaxi) 00469 END IF 00470 rij(i) = rij(i) + ddr 00471 dr = SQRT(SUM(rij(:)**2)) 00472 dfblock = -dfblock/ddr 00473 dsblock = -dsblock/ddr 00474 IF ( irow == iatom ) THEN 00475 foab = -SUM(dfblock*pblock) 00476 fow = SUM(dsblock*wblock) 00477 dsblock = -dsblock 00478 ELSE 00479 foab = SUM(dfblock*pblock) 00480 fow = -SUM(dsblock*wblock) 00481 END IF 00482 force_ab(i) = force_ab(i) + foab 00483 force_w(i) = force_w(i) + fow 00484 IF ( dftb_control%self_consistent ) THEN 00485 CPPrecondition(ASSOCIATED(dsblocks(i+1)%block),cp_failure_level,routineP,error,failure) 00486 dsblocks(i+1)%block = dsblocks(i+1)%block + dsblock 00487 END IF 00488 ENDDO 00489 IF ( use_virial ) THEN 00490 IF ( irow == iatom ) THEN 00491 f0 = -1._dp 00492 ELSE 00493 f0 = 1._dp 00494 END IF 00495 CALL virial_pair_force ( virial%pv_virial, -f0, force_ab, rij, error) 00496 CALL virial_pair_force ( virial%pv_virial, -f0, force_w, rij, error) 00497 IF ( qs_env%atprop%stress ) THEN 00498 f1 = 0.5_dp*f0 00499 CALL virial_pair_force (qs_env%atprop%atstress(:,:,iatom),-f1,force_ab,rij,error) 00500 CALL virial_pair_force (qs_env%atprop%atstress(:,:,iatom),-f1,force_w,rij,error) 00501 CALL virial_pair_force (qs_env%atprop%atstress(:,:,jatom),-f1,force_ab,rij,error) 00502 CALL virial_pair_force (qs_env%atprop%atstress(:,:,jatom),-f1,force_w,rij,error) 00503 END IF 00504 END IF 00505 DEALLOCATE (dfblock,dsblock,STAT=istat) 00506 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00507 END IF 00508 END IF 00509 00510 IF(calculate_forces .AND. (iatom/=jatom .OR. dr > 0.001_dp)) THEN 00511 atom_a = atom_of_kind(iatom) 00512 atom_b = atom_of_kind(jatom) 00513 IF ( irow == iatom ) force_ab = -force_ab 00514 IF ( irow == iatom ) force_w = -force_w 00515 force(ikind)%all_potential(:,atom_a) =& 00516 force(ikind)%all_potential(:,atom_a) - force_ab(:) 00517 force(jkind)%all_potential(:,atom_b) =& 00518 force(jkind)%all_potential(:,atom_b) + force_ab(:) 00519 force(ikind)%overlap(:,atom_a) =& 00520 force(ikind)%overlap(:,atom_a) - force_w(:) 00521 force(jkind)%overlap(:,atom_b) =& 00522 force(jkind)%overlap(:,atom_b) + force_w(:) 00523 END IF 00524 00525 ! gamma matrix 00526 IF ( dftb_control%self_consistent ) THEN 00527 hb_sr_damp = dftb_control%hb_sr_damp 00528 IF (hb_sr_damp) THEN 00529 ! short range correction enabled only when iatom XOR jatom are hydrogens 00530 hb_sr_damp = is_hydrogen(particle_set(iatom)%atomic_kind).NEQV.& 00531 is_hydrogen(particle_set(jatom)%atomic_kind) 00532 END IF 00533 ga = eta_a(0) 00534 gb = eta_b(0) 00535 gblock(1,1)= gblock(1,1) + gamma_rab_sr(dr,ga,gb,hb_sr_damp) 00536 IF(calculate_forces .AND. (iatom/=jatom .OR. dr > 0.001_dp)) THEN 00537 drp = dr + ddr 00538 drm = dr - ddr 00539 dgam = 0.5_dp*(gamma_rab_sr(drp,ga,gb,hb_sr_damp)-gamma_rab_sr(drm,ga,gb,hb_sr_damp))/ddr 00540 DO i=1,3 00541 CPPrecondition(ASSOCIATED(dgblocks(i+1)%block),cp_failure_level,routineP,error,failure) 00542 IF ( irow == iatom ) THEN 00543 dgblocks(i+1)%block(1,1)= dgblocks(i+1)%block(1,1) + dgam*rij(i)/dr 00544 ELSE 00545 dgblocks(i+1)%block(1,1)= dgblocks(i+1)%block(1,1) - dgam*rij(i)/dr 00546 END IF 00547 END DO 00548 END IF 00549 END IF 00550 00551 END IF 00552 00553 ! repulsive potential 00554 IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp) THEN 00555 erepij = 0._dp 00556 CALL urep_egr(rij,dr,erepij,force_rr,& 00557 n_urpoly,urep,spdim,s_cut,srep,spxr,scoeff,surr,calculate_forces) 00558 erep = erep + erepij 00559 IF(qs_env%atprop%energy) THEN 00560 qs_env%atprop%atecc(iatom) = qs_env%atprop%atecc(iatom) + 0.5_dp*erepij 00561 qs_env%atprop%atecc(jatom) = qs_env%atprop%atecc(jatom) + 0.5_dp*erepij 00562 END IF 00563 IF(calculate_forces .AND. (iatom/=jatom .OR. dr > 0.001_dp)) THEN 00564 atom_a = atom_of_kind(iatom) 00565 atom_b = atom_of_kind(jatom) 00566 force(ikind)%repulsive(:,atom_a) =& 00567 force(ikind)%repulsive(:,atom_a) - force_rr(:) 00568 force(jkind)%repulsive(:,atom_b) =& 00569 force(jkind)%repulsive(:,atom_b) + force_rr(:) 00570 IF ( use_virial ) THEN 00571 CALL virial_pair_force ( virial%pv_virial, -1._dp, force_rr, rij, error) 00572 IF(qs_env%atprop%stress) THEN 00573 CALL virial_pair_force(qs_env%atprop%atstress(:,:,iatom),-0.5_dp,force_rr,rij,error) 00574 CALL virial_pair_force(qs_env%atprop%atstress(:,:,jatom),-0.5_dp,force_rr,rij,error) 00575 END IF 00576 END IF 00577 END IF 00578 END IF 00579 00580 END DO 00581 CALL neighbor_list_iterator_release(nl_iterator) 00582 00583 IF ( dftb_control%self_consistent ) THEN 00584 DO i=1,SIZE(gamma_matrix) 00585 CALL cp_dbcsr_finalize(gamma_matrix(i)%matrix,error=error) 00586 ENDDO 00587 CALL set_qs_env(qs_env=qs_env,gamma_matrix=gamma_matrix,error=error) 00588 ENDIF 00589 DO i=1,SIZE(matrix_s) 00590 CALL cp_dbcsr_finalize(matrix_s(i)%matrix,error=error) 00591 ENDDO 00592 DO i=1,SIZE(matrix_h) 00593 CALL cp_dbcsr_finalize(matrix_h(i)%matrix,error=error) 00594 ENDDO 00595 00596 ! set repulsive energy 00597 CALL mp_sum(erep,para_env%group) 00598 energy%repulsive = erep 00599 00600 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00601 qs_env%input,"DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN",error=error),cp_p_file)) THEN 00602 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN",& 00603 extension=".Log",error=error) 00604 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1)%matrix,4,6,qs_env,para_env,& 00605 output_unit=iw,error=error) 00606 00607 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00608 qs_env%input,"DFT%PRINT%AO_MATRICES/DERIVATIVES",error=error),cp_p_file)) THEN 00609 DO i=2,SIZE(matrix_h) 00610 CALL cp_dbcsr_write_sparse_matrix(matrix_h(i)%matrix,4,6,qs_env,para_env,& 00611 output_unit=iw,error=error) 00612 END DO 00613 END IF 00614 00615 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 00616 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", error=error) 00617 END IF 00618 00619 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00620 qs_env%input,"DFT%PRINT%AO_MATRICES/OVERLAP",error=error),cp_p_file)) THEN 00621 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/OVERLAP",& 00622 extension=".Log",error=error) 00623 CALL cp_dbcsr_write_sparse_matrix(matrix_s(1)%matrix,4,6,qs_env,para_env,& 00624 output_unit=iw,error=error) 00625 00626 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00627 qs_env%input,"DFT%PRINT%AO_MATRICES/DERIVATIVES",error=error),cp_p_file)) THEN 00628 DO i=2,SIZE(matrix_s) 00629 CALL cp_dbcsr_write_sparse_matrix(matrix_s(i)%matrix,4,6,qs_env,para_env,& 00630 output_unit=iw,error=error) 00631 END DO 00632 END IF 00633 00634 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 00635 "DFT%PRINT%AO_MATRICES/OVERLAP", error=error) 00636 END IF 00637 00638 IF(calculate_forces) THEN 00639 IF (SIZE(matrix_p) == 2) THEN 00640 CALL cp_dbcsr_add(matrix_p(1)%matrix,matrix_p(2)%matrix,alpha_scalar=1.0_dp,& 00641 beta_scalar=-1.0_dp,error=error) 00642 END IF 00643 DEALLOCATE(atom_of_kind,STAT=istat) 00644 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00645 END IF 00646 00647 CALL timestop(handle) 00648 00649 END SUBROUTINE build_dftb_matrices 00650 00651 ! ***************************************************************************** 00652 SUBROUTINE build_dftb_ks_matrix (ks_env,qs_env,ks_matrix,rho,energy,& 00653 calculate_forces,just_energy,error) 00654 00655 TYPE(qs_ks_env_type), POINTER :: ks_env 00656 TYPE(qs_environment_type), POINTER :: qs_env 00657 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00658 POINTER :: ks_matrix 00659 TYPE(qs_rho_type), POINTER :: rho 00660 TYPE(qs_energy_type), POINTER :: energy 00661 LOGICAL, INTENT(in) :: calculate_forces, just_energy 00662 TYPE(cp_error_type), INTENT(inout) :: error 00663 00664 CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb_ks_matrix', 00665 routineP = moduleN//':'//routineN 00666 00667 INTEGER :: atom_a, handle, iatom, ikind, 00668 ispin, istat, natom, nkind, 00669 nspins, output_unit 00670 LOGICAL :: failure 00671 REAL(KIND=dp) :: zeff 00672 REAL(KIND=dp), DIMENSION(:), POINTER :: mcharge, occupation_numbers 00673 REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges 00674 TYPE(atomic_kind_type), DIMENSION(:), 00675 POINTER :: atomic_kind_set 00676 TYPE(atomic_kind_type), POINTER :: atomic_kind 00677 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00678 POINTER :: matrix_h, matrix_p, matrix_s, 00679 mo_derivs 00680 TYPE(cp_dbcsr_type), POINTER :: mo_coeff 00681 TYPE(cp_logger_type), POINTER :: logger 00682 TYPE(cp_para_env_type), POINTER :: para_env 00683 TYPE(dft_control_type), POINTER :: dft_control 00684 TYPE(mo_set_p_type), DIMENSION(:), 00685 POINTER :: mo_array 00686 TYPE(particle_type), DIMENSION(:), 00687 POINTER :: particle_set 00688 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind 00689 TYPE(section_vals_type), POINTER :: scf_section 00690 00691 CALL timeset(routineN,handle) 00692 00693 NULLIFY(dft_control, logger, scf_section) 00694 NULLIFY(matrix_p, particle_set) 00695 00696 energy%hartree = 0._dp 00697 00698 logger => cp_error_get_logger(error) 00699 00700 failure=.FALSE. 00701 CPPrecondition(ASSOCIATED(ks_env),cp_failure_level,routineP,error,failure) 00702 CPPrecondition(ks_env%ref_count>0,cp_failure_level,routineP,error,failure) 00703 00704 IF ( .NOT. failure ) THEN 00705 00706 CALL get_qs_env(qs_env=qs_env,& 00707 dft_control=dft_control,& 00708 atomic_kind_set=atomic_kind_set,& 00709 matrix_h=matrix_h,& 00710 para_env=para_env,error=error) 00711 00712 scf_section => section_vals_get_subs_vals(qs_env%input,"DFT%SCF",error=error) 00713 nspins=dft_control%nspins 00714 CPPrecondition(ASSOCIATED(matrix_h),cp_failure_level,routineP,error,failure) 00715 CPPrecondition(ASSOCIATED(rho),cp_failure_level,routineP,error,failure) 00716 CPPrecondition(SIZE(ks_matrix)>0,cp_failure_level,routineP,error,failure) 00717 00718 DO ispin=1,nspins 00719 ! copy the core matrix into the fock matrix 00720 CALL cp_dbcsr_copy(ks_matrix(ispin)%matrix,matrix_h(1)%matrix,error=error) 00721 END DO 00722 00723 IF ( dft_control%qs_control%dftb_control%self_consistent ) THEN 00724 ! Mulliken charges 00725 CALL get_qs_env(qs_env=qs_env,& 00726 particle_set=particle_set,& 00727 matrix_s=matrix_s,error=error) 00728 matrix_p => rho%rho_ao 00729 natom=SIZE(particle_set) 00730 ALLOCATE(charges(natom,nspins),STAT=istat) 00731 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00732 ! 00733 CALL mulliken_charges(matrix_p,matrix_s(1)%matrix,para_env,charges,error=error) 00734 ! 00735 ALLOCATE(mcharge(natom),STAT=istat) 00736 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00737 nkind = SIZE(atomic_kind_set) 00738 DO ikind=1,nkind 00739 atomic_kind => atomic_kind_set(ikind) 00740 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00741 natom=natom,dftb_parameter=dftb_kind) 00742 CALL get_dftb_atom_param(dftb_kind,zeff=zeff) 00743 DO iatom=1,natom 00744 atom_a = atomic_kind%atom_list(iatom) 00745 mcharge(atom_a) = zeff - SUM(charges(atom_a,1:nspins)) 00746 END DO 00747 END DO 00748 DEALLOCATE(charges,STAT=istat) 00749 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00750 00751 CALL build_dftb_coulomb(qs_env,ks_matrix,rho,mcharge,energy,& 00752 calculate_forces,just_energy,error) 00753 00754 DEALLOCATE(mcharge,STAT=istat) 00755 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00756 00757 END IF 00758 00759 energy%total = energy%core + energy%hartree + energy%qmmm_el + & 00760 energy%repulsive + energy%dispersion 00761 00762 output_unit=cp_print_key_unit_nr(logger,scf_section,"PRINT%DETAILED_ENERGY",& 00763 extension=".scfLog",error=error) 00764 IF (output_unit>0) THEN 00765 WRITE (UNIT=output_unit,FMT="(/,(T9,A,T60,F20.10))")& 00766 "Repulsive pair potential energy: ",energy%repulsive,& 00767 "Zeroth order Hamiltonian energy: ",energy%core,& 00768 "Charge fluctuation energy: ",energy%hartree,& 00769 "London dispersion energy: ",energy%dispersion 00770 IF (qs_env%qmmm) THEN 00771 WRITE (UNIT=output_unit,FMT="(T3,A,T60,F20.10)")& 00772 "QM/MM Electrostatic energy: ",energy%qmmm_el 00773 END IF 00774 END IF 00775 CALL cp_print_key_finished_output(output_unit,logger,scf_section,& 00776 "PRINT%DETAILED_ENERGY", error=error) 00777 ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C (plus occupation numbers) 00778 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN 00779 CALL get_qs_env(qs_env,mo_derivs=mo_derivs,mos=mo_array,error=error) 00780 DO ispin=1,SIZE(mo_derivs) 00781 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set,& 00782 mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers ) 00783 IF(.NOT.mo_array(ispin)%mo_set%use_mo_coeff_b) THEN 00784 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 00785 ENDIF 00786 CALL cp_dbcsr_multiply('n','n',1.0_dp,ks_matrix(ispin)%matrix,mo_coeff,& 00787 0.0_dp,mo_derivs(ispin)%matrix, error=error) 00788 ENDDO 00789 ENDIF 00790 00791 END IF 00792 00793 CALL timestop(handle) 00794 00795 END SUBROUTINE build_dftb_ks_matrix 00796 00797 ! ***************************************************************************** 00798 SUBROUTINE setup_matrices(qs_env,calculate_forces,error) 00799 00800 TYPE(qs_environment_type), POINTER :: qs_env 00801 LOGICAL :: calculate_forces 00802 TYPE(cp_error_type), INTENT(inout) :: error 00803 00804 CHARACTER(len=*), PARAMETER :: routineN = 'setup_matrices', 00805 routineP = moduleN//':'//routineN 00806 00807 INTEGER :: i, istat, natom, 00808 neighbor_list_id, nkind, 00809 nmat, nsgf 00810 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 00811 INTEGER, DIMENSION(:), POINTER :: rbs 00812 LOGICAL :: failure 00813 TYPE(array_i1d_obj) :: row_blk_sizes 00814 TYPE(atomic_kind_type), DIMENSION(:), 00815 POINTER :: atomic_kind_set 00816 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00817 POINTER :: matrix_h, matrix_s 00818 TYPE(dbcsr_distribution_obj), POINTER :: dbcsr_dist 00819 TYPE(dft_control_type), POINTER :: dft_control 00820 TYPE(dftb_control_type), POINTER :: dftb_control 00821 TYPE(neighbor_list_set_p_type), 00822 DIMENSION(:), POINTER :: sab_orb 00823 TYPE(particle_type), DIMENSION(:), 00824 POINTER :: particle_set 00825 00826 NULLIFY(matrix_s, matrix_h, particle_set, sab_orb, atomic_kind_set) 00827 00828 CALL get_qs_env(qs_env=qs_env,& 00829 matrix_s=matrix_s,& 00830 matrix_h=matrix_h,& 00831 atomic_kind_set=atomic_kind_set,& 00832 dft_control=dft_control,& 00833 particle_set=particle_set,& 00834 sab_orb=sab_orb,& 00835 dbcsr_dist=dbcsr_dist,& 00836 neighbor_list_id=neighbor_list_id, & 00837 error=error) 00838 00839 dftb_control => dft_control%qs_control%dftb_control 00840 00841 nkind = SIZE(atomic_kind_set) 00842 natom = SIZE(particle_set) 00843 00844 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,nsgf=nsgf) 00845 00846 ALLOCATE (first_sgf(natom),STAT=istat) 00847 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00848 ALLOCATE (last_sgf(natom),STAT=istat) 00849 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00850 00851 CALL get_particle_set(particle_set=particle_set,& 00852 first_sgf=first_sgf,& 00853 last_sgf=last_sgf,error=error) 00854 00855 IF ( dftb_control%self_consistent .AND. calculate_forces ) THEN 00856 ! we have to store the derivative overlap matrices 00857 nmat = 4 00858 ELSE 00859 nmat = 1 00860 END IF 00861 00862 ALLOCATE (rbs(natom), STAT=istat) 00863 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00864 CALL convert_offsets_to_sizes (first_sgf, rbs, last_sgf) 00865 CALL array_nullify (row_blk_sizes) 00866 CALL array_new (row_blk_sizes, rbs, gift=.TRUE.) 00867 00868 CALL cp_dbcsr_allocate_matrix_set(matrix_s,nmat,error=error) 00869 CALL cp_dbcsr_allocate_matrix_set(matrix_h,1,error=error) 00870 00871 ALLOCATE(matrix_s(1)%matrix) 00872 CALL cp_dbcsr_init(matrix_s(1)%matrix, error=error) 00873 CALL cp_dbcsr_create(matrix=matrix_s(1)%matrix, & 00874 name="OVERLAP MATRIX DFTB", & 00875 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric,& 00876 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00877 nblks=0, nze=0, mutable_work=.TRUE., & 00878 error=error) 00879 00880 DEALLOCATE (first_sgf,STAT=istat) 00881 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00882 DEALLOCATE (last_sgf,STAT=istat) 00883 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00884 00885 ALLOCATE(matrix_h(1)%matrix) 00886 CALL cp_dbcsr_init(matrix_h(1)%matrix, error=error) 00887 CALL cp_dbcsr_create(matrix=matrix_h(1)%matrix, & 00888 name="CORE HAMILTONIAN MATRIX DFTB", & 00889 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric,& 00890 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00891 nblks=0, nze=0, mutable_work=.TRUE., & 00892 error=error) 00893 00894 IF ( nmat > 1 ) THEN 00895 DO i=2,nmat 00896 ALLOCATE(matrix_s(i)%matrix) 00897 CALL cp_dbcsr_init(matrix_s(i)%matrix, error=error) 00898 CALL cp_dbcsr_create(matrix=matrix_s(i)%matrix, & 00899 name="OVERLAP DERIVATIVE MATRIX DFTB", & 00900 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric,& 00901 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00902 nblks=0, nze=0, mutable_work=.TRUE., & 00903 error=error) 00904 END DO 00905 END IF 00906 00907 CALL array_release (row_blk_sizes) 00908 00909 CALL set_qs_env(qs_env=qs_env,matrix_s=matrix_s,error=error) 00910 CALL set_qs_env(qs_env=qs_env,matrix_h=matrix_h,error=error) 00911 00912 END SUBROUTINE setup_matrices 00913 00914 ! ***************************************************************************** 00918 SUBROUTINE getskz(slakotab,skpar,dx,ngrd,ngrdcut,dgrd,llm) 00919 REAL(dp), INTENT(in) :: slakotab(:,:), dx 00920 INTEGER, INTENT(in) :: ngrd, ngrdcut 00921 REAL(dp), INTENT(in) :: dgrd 00922 INTEGER, INTENT(in) :: llm 00923 REAL(dp), INTENT(out) :: skpar(llm) 00924 00925 INTEGER :: clgp 00926 00927 skpar = 0._dp 00928 ! 00929 ! Determine closest grid point 00930 ! 00931 clgp = NINT(dx/dgrd) 00932 ! 00933 ! Screen elements which are too far away 00934 ! 00935 IF (clgp > ngrdcut) RETURN 00936 ! 00937 ! The grid point is either contained in the table --> matrix element 00938 ! can be interpolated, or it is outside the table --> matrix element 00939 ! needs to be extrapolated. 00940 ! 00941 IF (clgp > ngrd) THEN 00942 ! 00943 ! Extrapolate external matrix elements if table does not finish with zero 00944 ! 00945 CALL extrapol(slakotab,skpar,dx,ngrd,dgrd,llm) 00946 ELSE 00947 ! 00948 ! Interpolate tabulated matrix elements 00949 ! 00950 CALL interpol(slakotab,skpar,dx,ngrd,dgrd,llm,clgp) 00951 END IF 00952 END SUBROUTINE getskz 00953 00954 ! ***************************************************************************** 00955 SUBROUTINE interpol(slakotab,skpar,dx,ngrd,dgrd,llm,clgp) 00956 REAL(dp), INTENT(in) :: slakotab(:,:), dx 00957 INTEGER, INTENT(in) :: ngrd 00958 REAL(dp), INTENT(in) :: dgrd 00959 INTEGER, INTENT(in) :: llm 00960 REAL(dp), INTENT(out) :: skpar(llm) 00961 INTEGER, INTENT(in) :: clgp 00962 00963 INTEGER :: fgpm, k, l, lgpm 00964 REAL(dp) :: error, xa(max_inter), 00965 ya(max_inter) 00966 00967 lgpm = MIN(clgp+max_inter/2,ngrd) 00968 fgpm = lgpm - max_inter + 1 00969 DO k = 0,max_inter-1 00970 xa(k+1) = (fgpm+k)*dgrd 00971 END DO 00972 ! 00973 ! Interpolate matrix elements for all orbitals 00974 ! 00975 DO l = 1, llm 00976 ! 00977 ! Read SK parameters from table 00978 ! 00979 ya(1:max_inter) = slakotab(fgpm:lgpm,l) 00980 CALL polint(xa,ya,max_inter,dx,skpar(l),error) 00981 END DO 00982 END SUBROUTINE interpol 00983 00984 ! ***************************************************************************** 00985 SUBROUTINE extrapol(slakotab,skpar,dx,ngrd,dgrd,llm) 00986 REAL(dp), INTENT(in) :: slakotab(:,:), dx 00987 INTEGER, INTENT(in) :: ngrd 00988 REAL(dp), INTENT(in) :: dgrd 00989 INTEGER, INTENT(in) :: llm 00990 REAL(dp), INTENT(out) :: skpar(llm) 00991 00992 INTEGER :: fgp, k, l, lgp, ntable, nzero 00993 REAL(dp) :: error, xa(max_extra), 00994 ya(max_extra) 00995 00996 nzero = max_extra/3 00997 ntable = max_extra-nzero 00998 ! 00999 ! Get the three last distances from the table 01000 ! 01001 DO k = 1,ntable 01002 xa(k) = (ngrd-(max_extra-3)+k)*dgrd 01003 END DO 01004 DO k = 1,nzero 01005 xa(ntable+k) = (ngrd+k-1)*dgrd + slako_d0 01006 ya(ntable+k) = 0.0 01007 END DO 01008 ! 01009 ! Extrapolate matrix elements for all orbitals 01010 ! 01011 DO l = 1,llm 01012 ! 01013 ! Read SK parameters from table 01014 ! 01015 fgp = ngrd + 1 - (max_extra-3) 01016 lgp = ngrd 01017 ya(1:max_extra-3) = slakotab(fgp:lgp,l) 01018 CALL polint(xa,ya,max_extra,dx,skpar(l),error) 01019 END DO 01020 END SUBROUTINE extrapol 01021 01022 ! ***************************************************************************** 01038 SUBROUTINE turnsk(mat,skab1,skab2,dxv,dx,lmaxa,lmaxb) 01039 REAL(dp), INTENT(inout) :: mat(:,:) 01040 REAL(dp), INTENT(in) :: skab1(:), skab2(:), dxv(3), dx 01041 INTEGER, INTENT(in) :: lmaxa, lmaxb 01042 01043 INTEGER :: lmaxab, minlmaxab 01044 REAL(dp) :: rinv, rr(6), rr2(6) 01045 01046 lmaxab = MAX(lmaxa,lmaxb) 01047 ! Determine l quantum limits. 01048 IF (lmaxab.gt.2) STOP 'lmax=2' 01049 minlmaxab = MIN(lmaxa,lmaxb) 01050 ! 01051 ! s-s interaction 01052 ! 01053 CALL skss(skab1,mat) 01054 ! 01055 IF (lmaxab.le.0) RETURN 01056 ! 01057 rr2(1:3) = dxv(1:3)**2 01058 rr(1:3) = dxv(1:3) 01059 rinv = 1.0_dp/dx 01060 ! 01061 rr(1:3) = rr(1:3)*rinv 01062 rr(4:6) = rr(1:3) 01063 rr2(1:3) = rr2(1:3)*rinv**2 01064 rr2(4:6) = rr2(1:3) 01065 ! 01066 ! s-p, p-s and p-p interaction 01067 ! 01068 IF (minlmaxab.ge.1) THEN 01069 CALL skpp(skab1,mat,iptr(:,:,:,lmaxa,lmaxb)) 01070 CALL sksp(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01071 CALL sksp(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01072 ELSE 01073 IF (lmaxb.ge.1) THEN 01074 CALL sksp(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01075 ELSE 01076 CALL sksp(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01077 END IF 01078 END IF 01079 ! 01080 ! If there is only s-p interaction we have finished 01081 ! 01082 IF (lmaxab.le.1) RETURN 01083 ! 01084 ! at least one atom has d functions 01085 ! 01086 IF (minlmaxab.eq.2) THEN 01087 ! 01088 ! in case both atoms have d functions 01089 ! 01090 CALL skdd(skab1,mat,iptr(:,:,:,lmaxa,lmaxb)) 01091 CALL sksd(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01092 CALL sksd(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01093 CALL skpd(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01094 CALL skpd(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01095 ELSE 01096 ! 01097 ! One atom has d functions, the other has s or s and p functions 01098 ! 01099 IF (lmaxa.eq.0) THEN 01100 ! 01101 ! atom b has d, the atom a only s functions 01102 ! 01103 CALL sksd(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01104 ELSE IF (lmaxa.eq.1) THEN 01105 ! 01106 ! atom b has d, the atom a s and p functions 01107 ! 01108 CALL sksd(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01109 CALL skpd(skab2,mat,iptr(:,:,:,lmaxa,lmaxb),.TRUE.) 01110 ELSE 01111 ! 01112 ! atom a has d functions 01113 ! 01114 IF (lmaxb.eq.0) THEN 01115 ! 01116 ! atom a has d, atom b has only s functions 01117 ! 01118 CALL sksd(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01119 ELSE 01120 ! 01121 ! atom a has d, atom b has s and p functions 01122 ! 01123 CALL sksd(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01124 CALL skpd(skab1,mat,iptr(:,:,:,lmaxa,lmaxb),.FALSE.) 01125 END IF 01126 END IF 01127 END IF 01128 ! 01129 CONTAINS 01130 ! 01131 ! The subroutines to turn the matrix elements are taken as internal subroutines 01132 ! as it is beneficial to inline them. 01133 ! 01134 ! They are both turning the matrix elements and placing them appropriately 01135 ! into the matrix block 01136 ! 01137 ! ***************************************************************************** 01141 SUBROUTINE skss(skpar,mat) 01142 REAL(dp), INTENT(in) :: skpar(:) 01143 REAL(dp), INTENT(inout) :: mat(:,:) 01144 01145 mat(1,1) = mat(1,1) + skpar(1) 01146 ! 01147 END SUBROUTINE skss 01148 01149 ! ***************************************************************************** 01153 SUBROUTINE sksp(skpar,mat,ind,transposed) 01154 REAL(dp), INTENT(in) :: skpar(:) 01155 REAL(dp), INTENT(inout) :: mat(:,:) 01156 INTEGER, INTENT(in) :: ind(0:,0:,0:) 01157 LOGICAL, INTENT(in) :: transposed 01158 01159 INTEGER :: l 01160 REAL(dp) :: skp 01161 01162 skp = skpar(ind(1,0,0)) 01163 IF (transposed) THEN 01164 DO l = 1,3 01165 mat(1,l+1) = mat(1,l+1) + rr(l)*skp 01166 END DO 01167 ELSE 01168 DO l = 1,3 01169 mat(l+1,1) = mat(l+1,1) - rr(l)*skp 01170 END DO 01171 END IF 01172 ! 01173 END SUBROUTINE sksp 01174 01175 ! ***************************************************************************** 01176 SUBROUTINE skpp(skpar,mat,ind) 01177 REAL(dp), INTENT(in) :: skpar(:) 01178 REAL(dp), INTENT(inout) :: mat(:,:) 01179 INTEGER, INTENT(in) :: ind(0:,0:,0:) 01180 01181 INTEGER :: ii, ir, is, k, l 01182 REAL(dp) :: epp(6), matel(6), skppp, skpps 01183 01184 epp(1:3) = rr2(1:3) 01185 DO l = 1,3 01186 epp(l+3) = rr(l)*rr(l+1) 01187 END DO 01188 skppp = skpar(ind(1,1,1)) 01189 skpps = skpar(ind(1,1,0)) 01190 ! 01191 DO l = 1,3 01192 matel(l) = epp(l)*skpps + (1._dp-epp(l))*skppp 01193 END DO 01194 DO l = 4,6 01195 matel(l) = epp(l)*(skpps - skppp) 01196 END DO 01197 ! 01198 DO ir = 1,3 01199 DO is = 1,ir-1 01200 ii = ir - is 01201 k = 3*ii-(ii*(ii-1))/2+is 01202 mat(is+1,ir+1) = mat(is+1,ir+1) + matel(k) 01203 mat(ir+1,is+1) = mat(ir+1,is+1) + matel(k) 01204 END DO 01205 mat(ir+1,ir+1) = mat(ir+1,ir+1) + matel(ir) 01206 END DO 01207 END SUBROUTINE skpp 01208 01209 ! ***************************************************************************** 01210 SUBROUTINE sksd(skpar,mat,ind,transposed) 01211 REAL(dp), INTENT(in) :: skpar(:) 01212 REAL(dp), INTENT(inout) :: mat(:,:) 01213 INTEGER, INTENT(in) :: ind(0:,0:,0:) 01214 LOGICAL, INTENT(in) :: transposed 01215 01216 INTEGER :: l 01217 REAL(dp) :: d4, d5, es(5), r3, sksds 01218 01219 sksds = skpar(ind(2,0,0)) 01220 r3 = SQRT(3._dp) 01221 d4 = rr2(3) - 0.5_dp*(rr2(1)+rr2(2)) 01222 d5 = rr2(1) - rr2(2) 01223 ! 01224 DO l = 1,3 01225 es(l) = r3*rr(l)*rr(l+1) 01226 END DO 01227 es(4) = 0.5_dp*r3*d5 01228 es(5) = d4 01229 ! 01230 IF (transposed) THEN 01231 DO l = 1,5 01232 mat(1,l+4) = mat(1,l+4) + es(l)*sksds 01233 END DO 01234 ELSE 01235 DO l = 1,5 01236 mat(l+4,1) = mat(l+4,1) + es(l)*sksds 01237 END DO 01238 END IF 01239 END SUBROUTINE sksd 01240 01241 ! ***************************************************************************** 01242 SUBROUTINE skpd(skpar,mat,ind,transposed) 01243 REAL(dp), INTENT(in) :: skpar(:) 01244 REAL(dp), INTENT(inout) :: mat(:,:) 01245 INTEGER, INTENT(in) :: ind(0:,0:,0:) 01246 LOGICAL, INTENT(in) :: transposed 01247 01248 INTEGER :: ir, is, k, l, m 01249 REAL(dp) :: d3, d4, d5, d6, dm(15), 01250 epd(13,2), r3, sktmp 01251 01252 r3 = SQRT(3.0_dp) 01253 d3 = rr2(1) + rr2(2) 01254 d4 = rr2(3) - 0.5_dp*d3 01255 d5 = rr2(1) - rr2(2) 01256 d6 = rr(1)*rr(2)*rr(3) 01257 DO l = 1,3 01258 epd(l,1) = r3*rr2(l)*rr(l+1) 01259 epd(l,2) = rr(l+1)*(1.0_dp-2._dp*rr2(l)) 01260 epd(l+4,1) = r3*rr2(l)*rr(l+2) 01261 epd(l+4,2) = rr(l+2)*(1.0_dp-2*rr2(l)) 01262 epd(l+7,1) = 0.5_dp*r3*rr(l)*d5 01263 epd(l+10,1) = rr(l)*d4 01264 END DO 01265 ! 01266 epd(4,1) = r3*d6 01267 epd(4,2) = -2._dp*d6 01268 epd(8,2) = rr(1)*(1.0_dp-d5) 01269 epd(9,2) = -rr(2)*(1.0_dp+d5) 01270 epd(10,2) = -rr(3)*d5 01271 epd(11,2) = -r3*rr(1)*rr2(3) 01272 epd(12,2) = -r3*rr(2)*rr2(3) 01273 epd(13,2) = r3*rr(3)*d3 01274 ! 01275 dm(1:15) = 0.0_dp 01276 ! 01277 DO m = 1,2 01278 sktmp = skpar(ind(2,1,m-1)) 01279 dm(1)=dm(1)+epd(1,m)*sktmp 01280 dm(2)=dm(2)+epd(6,m)*sktmp 01281 dm(3)=dm(3)+epd(4,m)*sktmp 01282 dm(5)=dm(5)+epd(2,m)*sktmp 01283 dm(6)=dm(6)+epd(7,m)*sktmp 01284 dm(7)=dm(7)+epd(5,m)*sktmp 01285 dm(9)=dm(9)+epd(3,m)*sktmp 01286 DO l = 8,13 01287 dm(l+2) = dm(l+2)+epd(l,m)*sktmp 01288 END DO 01289 END DO 01290 ! 01291 dm(4) = dm(3) 01292 dm(8) = dm(3) 01293 ! 01294 IF (transposed) THEN 01295 DO ir = 1,5 01296 DO is = 1,3 01297 k=3*(ir-1)+is 01298 mat(is+1,ir+4) = mat(is+1,ir+4) + dm(k) 01299 END DO 01300 END DO 01301 ELSE 01302 DO ir = 1,5 01303 DO is = 1,3 01304 k=3*(ir-1)+is 01305 mat(ir+4,is+1) = mat(ir+4,is+1) - dm(k) 01306 END DO 01307 END DO 01308 END IF 01309 ! 01310 END SUBROUTINE skpd 01311 01312 ! ***************************************************************************** 01313 SUBROUTINE skdd(skpar,mat,ind) 01314 REAL(dp), INTENT(in) :: skpar(:) 01315 REAL(dp), INTENT(inout) :: mat(:,:) 01316 INTEGER, INTENT(in) :: ind(0:,0:,0:) 01317 01318 INTEGER :: ii, ir, is, k, l, m 01319 REAL(dp) :: d3, d4, d5, dd(3), dm(15), 01320 e(15,3), r3 01321 01322 r3 = SQRT(3._dp) 01323 d3 = rr2(1) + rr2(2) 01324 d4 = rr2(3) - 0.5_dp*d3 01325 d5 = rr2(1) - rr2(2) 01326 DO l = 1,3 01327 e(l,1) = rr2(l)*rr2(l+1) 01328 e(l,2) = rr2(l) + rr2(l+1) - 4._dp*e(l,1) 01329 e(l,3) = rr2(l+2) + e(l,1) 01330 e(l,1) = 3._dp*e(l,1) 01331 END DO 01332 e(4,1) = d5**2 01333 e(4,2) = d3 - e(4,1) 01334 e(4,3) = rr2(3) + 0.25_dp*e(4,1) 01335 e(4,1) = 0.75_dp*e(4,1) 01336 e(5,1) = d4**2 01337 e(5,2) = 3._dp*rr2(3)*d3 01338 e(5,3) = 0.75_dp*d3**2 01339 dd(1) = rr(1)*rr(3) 01340 dd(2) = rr(2)*rr(1) 01341 dd(3) = rr(3)*rr(2) 01342 DO l = 1,2 01343 e(l+5,1) = 3._dp*rr2(l+1)*dd(l) 01344 e(l+5,2) = dd(l)*(1._dp-4._dp*rr2(l+1)) 01345 e(l+5,3) = dd(l)*(rr2(l+1)-1._dp) 01346 END DO 01347 e(8,1) = dd(1)*d5*1.5_dp 01348 e(8,2) = dd(1)*(1.0_dp-2.0_dp*d5) 01349 e(8,3) = dd(1)*(0.5_dp*d5-1.0_dp) 01350 e(9,1) = d5*0.5_dp*d4*r3 01351 e(9,2) = -d5*rr2(3)*r3 01352 e(9,3) = d5*0.25_dp*(1.0_dp+rr2(3))*r3 01353 e(10,1) = rr2(1)*dd(3)*3.0_dp 01354 e(10,2) = (0.25_dp-rr2(1))*dd(3)*4.0_dp 01355 e(10,3) = dd(3)*(rr2(1)-1.0_dp) 01356 e(11,1) = 1.5_dp*dd(3)*d5 01357 e(11,2) = -dd(3)*(1.0_dp+2.0_dp*d5) 01358 e(11,3) = dd(3)*(1.0_dp+0.5_dp*d5) 01359 e(13,3) = 0.5_dp*d5*dd(2) 01360 e(13,2) = -2.0_dp*dd(2)*d5 01361 e(13,1) = e(13,3)*3.0_dp 01362 e(12,1) = d4*dd(1)*r3 01363 e(14,1) = d4*dd(3)*r3 01364 e(15,1) = d4*dd(2)*r3 01365 e(15,2) = -2.0_dp*r3*dd(2)*rr2(3) 01366 e(15,3) = 0.5_dp*r3*(1.0_dp+rr2(3))*dd(2) 01367 e(14,2) = r3*dd(3)*(d3-rr2(3)) 01368 e(14,3) = -r3*0.5_dp*dd(3)*d3 01369 e(12,2) = r3*dd(1)*(d3-rr2(3)) 01370 e(12,3) = -r3*0.5_dp*dd(1)*d3 01371 ! 01372 dm(1:15) = 0._dp 01373 DO l = 1,15 01374 DO m = 1,3 01375 dm(l) = dm(l)+e(l,m)*skpar(ind(2,2,m-1)) 01376 END DO 01377 END DO 01378 ! 01379 DO ir = 1,5 01380 DO is = 1,ir-1 01381 ii = ir-is 01382 k = 5*ii-(ii*(ii-1))/2+is 01383 mat(ir+4,is+4) = mat(ir+4,is+4) + dm(k) 01384 mat(is+4,ir+4) = mat(is+4,ir+4) + dm(k) 01385 END DO 01386 mat(ir+4,ir+4) = mat(ir+4,ir+4) + dm(ir) 01387 END DO 01388 END SUBROUTINE skdd 01389 ! 01390 END SUBROUTINE turnsk 01391 01392 ! ***************************************************************************** 01393 SUBROUTINE polint(xa,ya,n,x,y,dy) 01394 INTEGER, INTENT(in) :: n 01395 REAL(dp), INTENT(in) :: ya(n), xa(n), x 01396 REAL(dp), INTENT(out) :: y, dy 01397 01398 INTEGER :: i, m, ns 01399 REAL(dp) :: c(n), d(n), den, dif, dift, 01400 ho, hp, w 01401 01402 ! 01403 ! 01404 01405 ns=1 01406 01407 dif=ABS(x-xa(1)) 01408 DO i = 1,n 01409 dift=ABS(x-xa(i)) 01410 IF (dift.lt.dif) THEN 01411 ns=i 01412 dif=dift 01413 ENDIF 01414 c(i)=ya(i) 01415 d(i)=ya(i) 01416 END DO 01417 ! 01418 y=ya(ns) 01419 ns=ns-1 01420 DO m = 1,n-1 01421 DO i = 1,n-m 01422 ho=xa(i)-x 01423 hp=xa(i+m)-x 01424 w=c(i+1)-d(i) 01425 den=ho-hp 01426 IF(den.eq.0.) STOP 'POLINT' 01427 den=w/den 01428 d(i)=hp*den 01429 c(i)=ho*den 01430 END DO 01431 IF (2*ns.lt.n-m)THEN 01432 dy=c(ns+1) 01433 ELSE 01434 dy=d(ns) 01435 ns=ns-1 01436 ENDIF 01437 y=y+dy 01438 END DO 01439 ! 01440 RETURN 01441 END SUBROUTINE polint 01442 01443 ! ***************************************************************************** 01444 SUBROUTINE urep_egr(rv,r,erep,derep,& 01445 n_urpoly,urep,spdim,s_cut,srep,spxr,scoeff,surr,dograd) 01446 01447 REAL(dp), INTENT(in) :: rv(3), r 01448 REAL(dp), INTENT(inout) :: erep, derep(3) 01449 INTEGER, INTENT(in) :: n_urpoly 01450 REAL(dp), INTENT(in) :: urep(:) 01451 INTEGER, INTENT(in) :: spdim 01452 REAL(dp), INTENT(in) :: s_cut, srep(3) 01453 REAL(dp), POINTER :: spxr(:,:), scoeff(:,:) 01454 REAL(dp), INTENT(in) :: surr(2) 01455 LOGICAL, INTENT(in) :: dograd 01456 01457 INTEGER :: ic, isp, jsp, nsp 01458 REAL(dp) :: de_z, rz 01459 01460 derep=0._dp 01461 de_z = 0._dp 01462 IF (n_urpoly > 0) THEN 01463 ! 01464 ! polynomial part 01465 ! 01466 rz = urep(1) - r 01467 IF (rz <= rtiny) RETURN 01468 DO ic = 2,n_urpoly 01469 erep = erep + urep(ic) * rz**(ic) 01470 END DO 01471 IF (dograd) THEN 01472 DO ic = 2,n_urpoly 01473 de_z = de_z - ic*urep(ic) * rz**(ic-1) 01474 END DO 01475 END IF 01476 ELSE IF (spdim > 0) THEN 01477 ! 01478 ! spline part 01479 ! 01480 ! This part is kind of proprietary Paderborn code and I won't reverse-engeneer 01481 ! everything in detail. What is obvious is documented. 01482 ! 01483 ! This part has 4 regions: 01484 ! a) very long range is screened 01485 ! b) short-range is extrapolated with e-functions 01486 ! ca) normal range is approximated with a spline 01487 ! cb) longer range is extrapolated with an higher degree spline 01488 ! 01489 IF (r > s_cut) RETURN ! screening (condition a) 01490 ! 01491 IF (r < spxr(1,1)) THEN 01492 ! a) short range 01493 erep = erep + EXP(-srep(1)*r + srep(2)) + srep(3) 01494 IF (dograd) de_z = de_z -srep(1)*EXP(-srep(1)*r + srep(2)) 01495 ELSE 01496 ! 01497 ! condition c). First determine between which places the spline is located: 01498 ! 01499 ispg: DO isp = 1,spdim ! condition ca) 01500 IF (r < spxr(isp,1)) CYCLE ispg ! distance is smaller than this spline range 01501 IF (r >= spxr(isp,2)) CYCLE ispg ! distance is larger than this spline range 01502 ! at this point we have found the correct spline interval 01503 rz = r - spxr(isp,1) 01504 IF (isp /= spdim) THEN 01505 nsp = 3 ! condition ca 01506 DO jsp = 0,nsp 01507 erep = erep + scoeff(isp,jsp+1)*rz**(jsp) 01508 END DO 01509 IF (dograd) THEN 01510 DO jsp = 1,nsp 01511 de_z = de_z + jsp*scoeff(isp,jsp+1)*rz**(jsp-1) 01512 END DO 01513 END IF 01514 ELSE 01515 nsp = 5 ! condition cb 01516 DO jsp = 0,nsp 01517 IF( jsp <= 3 ) THEN 01518 erep = erep + scoeff(isp,jsp+1)*rz**(jsp) 01519 ELSE 01520 erep = erep + surr(jsp-3)*rz**(jsp) 01521 ENDIF 01522 END DO 01523 IF (dograd) THEN 01524 DO jsp = 1,nsp 01525 IF( jsp <= 3 ) THEN 01526 de_z = de_z + jsp*scoeff(isp,jsp+1)*rz**(jsp-1) 01527 ELSE 01528 de_z = de_z + jsp*surr(jsp-3)*rz**(jsp-1) 01529 ENDIF 01530 END DO 01531 END IF 01532 END IF 01533 EXIT ispg 01534 END DO ispg 01535 END IF 01536 END IF 01537 ! 01538 IF (dograd) THEN 01539 IF ( r > 1.e-12_dp ) derep(1:3) = (de_z/r)*rv(1:3) 01540 END IF 01541 01542 END SUBROUTINE urep_egr 01543 01544 ! ***************************************************************************** 01553 FUNCTION gamma_rab_sr(r,ga,gb,sr_damp) RESULT(gamma) 01554 REAL(dp), INTENT(in) :: r, ga, gb 01555 LOGICAL, INTENT(in) :: sr_damp 01556 REAL(dp) :: gamma 01557 01558 REAL(dp) :: a, b, fac, g_sum 01559 01560 gamma = 0.0_dp 01561 a = 3.2_dp*ga ! 3.2 = 16/5 in Eq. 18 and ff. 01562 b = 3.2_dp*gb 01563 g_sum = a + b 01564 IF (g_sum.lt.tol_gamma) RETURN ! hardness screening 01565 IF (r < rtiny) THEN ! This is for short distances but non-onsite terms 01566 ! This gives also correct diagonal elements (a=b, r=0) 01567 gamma = 0.5_dp*(a*b/g_sum + (a*b)**2/g_sum**3) 01568 RETURN 01569 END IF 01570 ! 01571 ! distinguish two cases: Gamma's are very close, e.g. for the same atom type, 01572 ! and Gamma's are different 01573 ! 01574 IF (ABS(a-b) < rtiny) THEN 01575 fac = 1.6_dp*r*a*b/g_sum*(1.0_dp + a*b/g_sum**2) 01576 gamma = -(48.0_dp + 33._dp*fac + (9.0_dp + fac)*fac**2)*EXP(-fac)/(48._dp*r) 01577 ELSE 01578 gamma = -EXP(-a*r)*(0.5_dp*a*b**4/(a**2-b**2)**2 - & 01579 (b**6 - 3._dp*a**2*b**4)/(r*(a**2-b**2)**3)) - & ! a-> b 01580 EXP(-b*r)*(0.5_dp*b*a**4/(b**2-a**2)**2 - & 01581 (a**6 - 3._dp*b**2*a**4)/(r*(b**2-a**2)**3)) ! b-> a 01582 END IF 01583 ! 01584 ! damping function for better short range hydrogen bonds. 01585 ! functional form from Hu H. et al., J. Phys. Chem. A 2007, 111, 5685-5691 01586 ! according to Elstner M, Theor. Chem. Acc. 2006, 116, 316-325, 01587 ! this should only be applied to a-b pairs involving hydrogen. 01588 IF (sr_damp) THEN 01589 gamma = gamma * EXP(-(0.5_dp*(ga+gb))**4 * r*r) 01590 END IF 01591 END FUNCTION gamma_rab_sr 01592 01593 END MODULE qs_dftb_matrices 01594
1.7.3