|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00011 MODULE qs_fermi_contact 00012 00013 USE ai_fermi_contact, ONLY: fermi_contact 00014 USE atomic_kind_types, ONLY: atomic_kind_type,& 00015 get_atomic_kind,& 00016 get_atomic_kind_set 00017 USE basis_set_types, ONLY: gto_basis_set_p_type,& 00018 gto_basis_set_type 00019 USE block_p_types, ONLY: block_p_type 00020 USE cell_types, ONLY: cell_type,& 00021 pbc 00022 USE cp_dbcsr_interface, ONLY: cp_dbcsr_get_block_p 00023 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 00024 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00025 USE cp_output_handling, ONLY: cp_p_file,& 00026 cp_print_key_finished_output,& 00027 cp_print_key_should_output,& 00028 cp_print_key_unit_nr 00029 USE cp_para_types, ONLY: cp_para_env_type 00030 USE kinds, ONLY: dp,& 00031 dp_size,& 00032 int_size 00033 USE orbital_pointers, ONLY: init_orbital_pointers,& 00034 ncoset 00035 USE particle_types, ONLY: particle_type 00036 USE qs_environment_types, ONLY: get_qs_env,& 00037 qs_environment_type 00038 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 00039 neighbor_list_iterate,& 00040 neighbor_list_iterator_create,& 00041 neighbor_list_iterator_p_type,& 00042 neighbor_list_iterator_release,& 00043 neighbor_list_set_p_type 00044 USE termination, ONLY: stop_memory 00045 USE timings, ONLY: timeset,& 00046 timestop 00047 #include "cp_common_uses.h" 00048 00049 IMPLICIT NONE 00050 00051 PRIVATE 00052 00053 ! *** Global parameters *** 00054 00055 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fermi_contact' 00056 00057 ! *** Public subroutines *** 00058 00059 PUBLIC :: build_fermi_contact_matrix 00060 00061 CONTAINS 00062 00063 ! ***************************************************************************** 00070 00071 SUBROUTINE build_fermi_contact_matrix(qs_env,matrix_fc,rc,error) 00072 00073 TYPE(qs_environment_type), POINTER :: qs_env 00074 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00075 POINTER :: matrix_fc 00076 REAL(dp), DIMENSION(3), INTENT(IN) :: rc 00077 TYPE(cp_error_type), INTENT(inout) :: error 00078 00079 CHARACTER(len=*), PARAMETER :: routineN = 'build_fermi_contact_matrix', 00080 routineP = moduleN//':'//routineN 00081 00082 INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, istat, iw, 00083 jatom, jkind, jset, last_jatom, ldai, ldfc, maxco, maxlgto, maxsgf, 00084 natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb 00085 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, 00086 lb_min, npgfa, npgfb, nsgfa, 00087 nsgfb 00088 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 00089 LOGICAL :: found, new_atom_b 00090 REAL(KIND=dp) :: dab, rab2 00091 REAL(KIND=dp), ALLOCATABLE, 00092 DIMENSION(:, :) :: fcab, work 00093 REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc 00094 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 00095 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, 00096 zeta, zetb 00097 TYPE(atomic_kind_type), DIMENSION(:), 00098 POINTER :: atomic_kind_set 00099 TYPE(atomic_kind_type), POINTER :: atomic_kind 00100 TYPE(block_p_type), ALLOCATABLE, 00101 DIMENSION(:) :: fcint 00102 TYPE(cell_type), POINTER :: cell 00103 TYPE(cp_logger_type), POINTER :: logger 00104 TYPE(cp_para_env_type), POINTER :: para_env 00105 TYPE(gto_basis_set_p_type), 00106 DIMENSION(:), POINTER :: basis_set_list 00107 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 00108 TYPE(neighbor_list_iterator_p_type), 00109 DIMENSION(:), POINTER :: nl_iterator 00110 TYPE(neighbor_list_set_p_type), 00111 DIMENSION(:), POINTER :: sab_orb 00112 TYPE(particle_type), DIMENSION(:), 00113 POINTER :: particle_set 00114 00115 CALL timeset(routineN,handle) 00116 00117 NULLIFY(cell,sab_orb,atomic_kind_set,particle_set,para_env) 00118 NULLIFY(logger) 00119 00120 logger => cp_error_get_logger(error) 00121 00122 CALL get_qs_env(qs_env=qs_env,& 00123 atomic_kind_set=atomic_kind_set,& 00124 particle_set=particle_set,& 00125 para_env=para_env,& 00126 sab_orb=sab_orb,& 00127 cell=cell,& 00128 error=error) 00129 00130 nkind = SIZE(atomic_kind_set) 00131 natom = SIZE(particle_set) 00132 00133 ! *** Allocate work storage *** 00134 00135 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00136 maxco=maxco,& 00137 maxlgto=maxlgto,& 00138 maxsgf=maxsgf) 00139 00140 ldai = ncoset(maxlgto) 00141 CALL init_orbital_pointers(ldai) 00142 00143 ldfc = maxco 00144 ALLOCATE (fcab(ldfc,ldfc),STAT=istat) 00145 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00146 "fcab",ldfc*ldfc*dp_size) 00147 fcab(:,:) = 0.0_dp 00148 00149 ALLOCATE (work(maxco,maxsgf),STAT=istat) 00150 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00151 "work",maxco*maxsgf*dp_size) 00152 work(:,:) = 0.0_dp 00153 00154 ALLOCATE (fcint(1),STAT=istat) 00155 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00156 "fcint",1*int_size) 00157 NULLIFY(fcint(1)%block) 00158 00159 ALLOCATE (basis_set_list(nkind),STAT=istat) 00160 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00161 "basis_set_list",nkind) 00162 DO ikind=1,nkind 00163 atomic_kind => atomic_kind_set(ikind) 00164 CALL get_atomic_kind(atomic_kind=atomic_kind,orb_basis_set=basis_set_a) 00165 IF (ASSOCIATED(basis_set_a)) THEN 00166 basis_set_list(ikind)%gto_basis_set => basis_set_a 00167 ELSE 00168 NULLIFY(basis_set_list(ikind)%gto_basis_set) 00169 END IF 00170 END DO 00171 CALL neighbor_list_iterator_create(nl_iterator,sab_orb) 00172 DO WHILE (neighbor_list_iterate(nl_iterator)==0) 00173 CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,inode=inode,& 00174 iatom=iatom,jatom=jatom,r=rab) 00175 basis_set_a => basis_set_list(ikind)%gto_basis_set 00176 IF (.NOT.ASSOCIATED(basis_set_a)) CYCLE 00177 basis_set_b => basis_set_list(jkind)%gto_basis_set 00178 IF (.NOT.ASSOCIATED(basis_set_b)) CYCLE 00179 ra = pbc(particle_set(iatom)%r,cell) 00180 ! basis ikind 00181 first_sgfa => basis_set_a%first_sgf 00182 la_max => basis_set_a%lmax 00183 la_min => basis_set_a%lmin 00184 npgfa => basis_set_a%npgf 00185 nseta = basis_set_a%nset 00186 nsgfa => basis_set_a%nsgf_set 00187 rpgfa => basis_set_a%pgf_radius 00188 set_radius_a => basis_set_a%set_radius 00189 sphi_a => basis_set_a%sphi 00190 zeta => basis_set_a%zet 00191 ! basis jkind 00192 first_sgfb => basis_set_b%first_sgf 00193 lb_max => basis_set_b%lmax 00194 lb_min => basis_set_b%lmin 00195 npgfb => basis_set_b%npgf 00196 nsetb = basis_set_b%nset 00197 nsgfb => basis_set_b%nsgf_set 00198 rpgfb => basis_set_b%pgf_radius 00199 set_radius_b => basis_set_b%set_radius 00200 sphi_b => basis_set_b%sphi 00201 zetb => basis_set_b%zet 00202 00203 IF(inode==1) last_jatom = 0 00204 00205 rb = rab + ra 00206 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 00207 dab = SQRT(rab2) 00208 rac = pbc(ra,rc,cell) 00209 rbc = rac - rab 00210 00211 IF (jatom /= last_jatom) THEN 00212 new_atom_b = .TRUE. 00213 last_jatom = jatom 00214 ELSE 00215 new_atom_b = .FALSE. 00216 END IF 00217 00218 IF (new_atom_b) THEN 00219 IF (iatom <= jatom) THEN 00220 irow = iatom 00221 icol = jatom 00222 ELSE 00223 irow = jatom 00224 icol = iatom 00225 END IF 00226 00227 NULLIFY (fcint(1)%block) 00228 CALL cp_dbcsr_get_block_p(matrix=matrix_fc(1)%matrix,& 00229 row=irow,col=icol,BLOCK=fcint(1)%block,found=found) 00230 ENDIF 00231 00232 DO iset=1,nseta 00233 00234 ncoa = npgfa(iset)*ncoset(la_max(iset)) 00235 sgfa = first_sgfa(1,iset) 00236 00237 DO jset=1,nsetb 00238 00239 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 00240 00241 ncob = npgfb(jset)*ncoset(lb_max(jset)) 00242 sgfb = first_sgfb(1,jset) 00243 00244 ! *** Calculate the primitive fermi contact integrals *** 00245 00246 CALL fermi_contact(la_max(iset),la_min(iset),npgfa(iset),& 00247 rpgfa(:,iset),zeta(:,iset),& 00248 lb_max(jset),lb_min(jset),npgfb(jset),& 00249 rpgfb(:,jset),zetb(:,jset),& 00250 rac,rbc,dab,fcab,SIZE(fcab,1)) 00251 00252 ! *** Contraction step *** 00253 00254 CALL dgemm("N","N",ncoa,nsgfb(jset),ncob,& 00255 1.0_dp,fcab(1,1),SIZE(fcab,1),& 00256 sphi_b(1,sgfb),SIZE(sphi_b,1),& 00257 0.0_dp,work(1,1),SIZE(work,1)) 00258 00259 IF (iatom <= jatom) THEN 00260 00261 CALL dgemm("T","N",nsgfa(iset),nsgfb(jset),ncoa,& 00262 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 00263 work(1,1),SIZE(work,1),& 00264 1.0_dp,fcint(1)%block(sgfa,sgfb),& 00265 SIZE(fcint(1)%block,1)) 00266 00267 ELSE 00268 00269 CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,& 00270 1.0_dp,work(1,1),SIZE(work,1),& 00271 sphi_a(1,sgfa),SIZE(sphi_a,1),& 00272 1.0_dp,fcint(1)%block(sgfb,sgfa),& 00273 SIZE(fcint(1)%block,1)) 00274 00275 ENDIF 00276 00277 ENDDO 00278 00279 ENDDO 00280 00281 END DO 00282 CALL neighbor_list_iterator_release(nl_iterator) 00283 00284 ! *** Release work storage *** 00285 DEALLOCATE (basis_set_list,STAT=istat) 00286 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00287 "basis_set_list") 00288 00289 DEALLOCATE (fcab,STAT=istat) 00290 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00291 "fcab") 00292 00293 DEALLOCATE (work,STAT=istat) 00294 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00295 "work") 00296 00297 NULLIFY (fcint(1)%block) 00298 DEALLOCATE (fcint,STAT=istat) 00299 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00300 "fcint") 00301 00302 ! *** Print the Fermi contact matrix, if requested *** 00303 00304 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 00305 qs_env%input,"DFT%PRINT%AO_MATRICES/FERMI_CONTACT",error=error),cp_p_file)) THEN 00306 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/FERMI_CONTACT",& 00307 extension=".Log",error=error) 00308 CALL cp_dbcsr_write_sparse_matrix(matrix_fc(1)%matrix,4,6,qs_env,para_env,output_unit=iw,error=error) 00309 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 00310 "DFT%PRINT%AO_MATRICES/FERMI_CONTACT", error=error) 00311 END IF 00312 00313 CALL timestop(handle) 00314 00315 END SUBROUTINE build_fermi_contact_matrix 00316 00317 ! ***************************************************************************** 00318 00319 END MODULE qs_fermi_contact 00320
1.7.3