|
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 ! ***************************************************************************** 00007 MODULE atom_energy 00008 USE atom_electronic_structure, ONLY: calculate_atom 00009 USE atom_fit, ONLY: atom_fit_density 00010 USE atom_operators, ONLY: atom_int_release,& 00011 atom_int_setup,& 00012 atom_ppint_release,& 00013 atom_ppint_setup,& 00014 atom_relint_release,& 00015 atom_relint_setup 00016 USE atom_output, ONLY: atom_print_basis,& 00017 atom_print_info,& 00018 atom_print_method,& 00019 atom_print_orbitals,& 00020 atom_print_potential,& 00021 atom_write_pseudo_param 00022 USE atom_types, ONLY: & 00023 GTH_PSEUDO, atom_basis_type, atom_gthpot_type, atom_integrals, & 00024 atom_optimization_type, atom_orbitals, atom_p_type, & 00025 atom_potential_type, atom_state, atom_type, create_atom_orbs, & 00026 create_atom_type, init_atom_basis, init_atom_potential, & 00027 read_atom_opt_section, release_atom_basis, release_atom_potential, & 00028 release_atom_type, set_atom 00029 USE atom_utils, ONLY: atom_consistent_method,& 00030 atom_core_density,& 00031 atom_density,& 00032 atom_local_potential,& 00033 atom_set_occupation,& 00034 get_maxl_occ,& 00035 get_maxn_occ 00036 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00037 cp_print_key_unit_nr 00038 USE f77_blas 00039 USE input_constants, ONLY: do_analytic 00040 USE input_section_types, ONLY: section_vals_get,& 00041 section_vals_get_subs_vals,& 00042 section_vals_type,& 00043 section_vals_val_get 00044 USE kinds, ONLY: default_string_length,& 00045 dp 00046 USE lapack 00047 USE mathconstants, ONLY: dfac,& 00048 gamma1,& 00049 pi 00050 USE periodic_table, ONLY: nelem,& 00051 ptable 00052 USE timings, ONLY: timeset,& 00053 timestop 00054 #include "cp_common_uses.h" 00055 00056 IMPLICIT NONE 00057 PRIVATE 00058 PUBLIC :: atom_energy_opt 00059 00060 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_energy' 00061 00062 ! ***************************************************************************** 00063 00064 CONTAINS 00065 00066 ! ***************************************************************************** 00067 00068 SUBROUTINE atom_energy_opt(atom_section,error) 00069 TYPE(section_vals_type), POINTER :: atom_section 00070 TYPE(cp_error_type), INTENT(INOUT) :: error 00071 00072 CHARACTER(len=*), PARAMETER :: routineN = 'atom_energy_opt', 00073 routineP = moduleN//':'//routineN 00074 00075 CHARACTER(LEN=2) :: elem 00076 CHARACTER(LEN=default_string_length), 00077 DIMENSION(:), POINTER :: tmpstringlist 00078 INTEGER :: do_eric, do_erie, handle, i, ierr, im, in, iw, k, maxl, mb, 00079 method, mo, n_meth, n_rep, nder, num_gto, reltyp, zcore, zval, zz 00080 INTEGER, DIMENSION(0:3) :: maxn 00081 INTEGER, DIMENSION(:), POINTER :: cn 00082 LOGICAL :: eri_c, eri_e, failure, 00083 had_ae, had_pp, pp_calc 00084 REAL(KIND=dp) :: delta 00085 REAL(KIND=dp), DIMENSION(0:3, 10) :: pocc 00086 TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis 00087 TYPE(atom_integrals), POINTER :: ae_int, pp_int 00088 TYPE(atom_optimization_type) :: optimization 00089 TYPE(atom_orbitals), POINTER :: orbitals 00090 TYPE(atom_p_type), DIMENSION(:, :), 00091 POINTER :: atom_info 00092 TYPE(atom_potential_type), POINTER :: ae_pot, p_pot 00093 TYPE(atom_state), POINTER :: state 00094 TYPE(cp_logger_type), POINTER :: logger 00095 TYPE(section_vals_type), POINTER :: basis_section, method_section, 00096 opt_section, potential_section, powell_section, xc_section 00097 00098 failure = .FALSE. 00099 00100 CALL timeset(routineN,handle) 00101 00102 ! What atom do we calculate 00103 CALL section_vals_val_get(atom_section,"ATOMIC_NUMBER", i_val=zval, error=error) 00104 CALL section_vals_val_get(atom_section,"ELEMENT", c_val=elem, error=error) 00105 zz = 0 00106 DO i=1,nelem 00107 IF ( ptable(i)%symbol == elem ) THEN 00108 zz = i 00109 EXIT 00110 END IF 00111 END DO 00112 IF ( zz /= 1 ) zval = zz 00113 00114 ! read and set up inofrmation on the basis sets 00115 ALLOCATE(ae_basis,pp_basis,STAT=ierr) 00116 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00117 basis_section => section_vals_get_subs_vals(atom_section,"AE_BASIS",error=error) 00118 NULLIFY(ae_basis%grid) 00119 CALL init_atom_basis(ae_basis,basis_section,zval,"AE",error) 00120 NULLIFY(pp_basis%grid) 00121 basis_section => section_vals_get_subs_vals(atom_section,"PP_BASIS",error=error) 00122 CALL init_atom_basis(pp_basis,basis_section,zval,"PP",error) 00123 00124 ! print general and basis set information 00125 logger => cp_error_get_logger(error) 00126 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%PROGRAM_BANNER",extension=".log",error=error) 00127 IF(iw > 0) CALL atom_print_info(zval,"Atomic Energy Calculation",iw,error) 00128 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%PROGRAM_BANNER",error=error) 00129 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%BASIS_SET",extension=".log",error=error) 00130 IF(iw > 0) THEN 00131 CALL atom_print_basis(ae_basis,iw," All Electron Basis",error) 00132 CALL atom_print_basis(pp_basis,iw," Pseudopotential Basis",error) 00133 END IF 00134 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%BASIS_SET",error=error) 00135 00136 ! read and setup information on the pseudopotential 00137 NULLIFY(potential_section) 00138 potential_section => section_vals_get_subs_vals(atom_section,"POTENTIAL",error=error) 00139 ALLOCATE(ae_pot,p_pot,STAT=ierr) 00140 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00141 CALL init_atom_potential(p_pot,potential_section,zval,error) 00142 CALL init_atom_potential(ae_pot,potential_section,-1,error) 00143 00144 ! if the ERI's are calculated analytically, we have to precalculate them 00145 eri_c = .FALSE. 00146 CALL section_vals_val_get(atom_section,"COULOMB_INTEGRALS", i_val=do_eric, error=error) 00147 IF(do_eric==do_analytic) eri_c = .TRUE. 00148 eri_e = .FALSE. 00149 CALL section_vals_val_get(atom_section,"EXCHANGE_INTEGRALS", i_val=do_erie, error=error) 00150 IF(do_erie==do_analytic) eri_e = .TRUE. 00151 00152 ! information on the states to be calculated 00153 CALL section_vals_val_get(atom_section,"MAX_ANGULAR_MOMENTUM", i_val=maxl, error=error) 00154 maxn=0 00155 CALL section_vals_val_get(atom_section,"CALCULATE_STATES", i_vals=cn, error=error) 00156 DO in = 1, MIN(SIZE(cn),4) 00157 maxn(in-1) = cn(in) 00158 END DO 00159 DO in = 0, 3 00160 maxn(in) = MIN(maxn(in),ae_basis%nbas(in)) 00161 maxn(in) = MIN(maxn(in),pp_basis%nbas(in)) 00162 END DO 00163 00164 ! read optimization section 00165 opt_section => section_vals_get_subs_vals(atom_section,"OPTIMIZATION",error=error) 00166 CALL read_atom_opt_section(optimization,opt_section,error) 00167 00168 had_ae = .FALSE. 00169 had_pp = .FALSE. 00170 00171 ! Check for the total number of electron configurations to be calculated 00172 CALL section_vals_val_get(atom_section,"ELECTRON_CONFIGURATION", n_rep_val=n_rep, error=error) 00173 ! Check for the total number of method types to be calculated 00174 method_section => section_vals_get_subs_vals(atom_section,"METHOD",error=error) 00175 CALL section_vals_get(method_section,n_repetition=n_meth,error=error) 00176 00177 ! integrals 00178 ALLOCATE(ae_int, pp_int,STAT=ierr) 00179 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00180 00181 ALLOCATE(atom_info(n_rep,n_meth),STAT=ierr) 00182 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00183 00184 DO in = 1, n_rep 00185 DO im = 1, n_meth 00186 00187 NULLIFY(atom_info(in,im)%atom) 00188 CALL create_atom_type(atom_info(in,im)%atom,error) 00189 00190 atom_info(in,im)%atom%optimization = optimization 00191 00192 atom_info(in,im)%atom%z = zval 00193 xc_section => section_vals_get_subs_vals(method_section,"XC",i_rep_section=im,error=error) 00194 atom_info(in,im)%atom%xc_section => xc_section 00195 00196 ALLOCATE(state,STAT=ierr) 00197 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00198 00199 ! get the electronic configuration 00200 CALL section_vals_val_get(atom_section,"ELECTRON_CONFIGURATION", i_rep_val=in,& 00201 c_vals=tmpstringlist, error=error) 00202 00203 ! set occupations 00204 CALL atom_set_occupation(tmpstringlist,state%occ,state%occupation,state%multiplicity,error) 00205 state%maxl_occ = get_maxl_occ(state%occ) 00206 state%maxn_occ = get_maxn_occ(state%occ) 00207 00208 ! set number of states to be calculated 00209 state%maxl_calc = MAX(maxl,state%maxl_occ) 00210 state%maxl_calc = MIN(3,state%maxl_calc) 00211 state%maxn_calc = 0 00212 DO k=0,state%maxl_calc 00213 state%maxn_calc(k) = MAX(maxn(k),state%maxn_occ(k)) 00214 END DO 00215 00216 ! is there a pseudo potential 00217 pp_calc = ANY(INDEX(tmpstringlist(1:),"CORE") /= 0) 00218 IF ( pp_calc ) THEN 00219 ! get and set the core occupations 00220 CALL section_vals_val_get(atom_section,"CORE", c_vals=tmpstringlist, error=error) 00221 CALL atom_set_occupation(tmpstringlist,state%core,pocc,error=error) 00222 zcore = zval - SUM(state%core) 00223 CALL set_atom(atom_info(in,im)%atom,zcore=zcore,pp_calc=.TRUE.,error=error) 00224 ELSE 00225 state%core=0._dp 00226 CALL set_atom(atom_info(in,im)%atom,zcore=zval,pp_calc=.FALSE.,error=error) 00227 END IF 00228 00229 CALL section_vals_val_get(method_section,"METHOD_TYPE",i_val=method,i_rep_section=im,error=error) 00230 CALL section_vals_val_get(method_section,"RELATIVISTIC",i_val=reltyp,i_rep_section=im,error=error) 00231 CALL set_atom(atom_info(in,im)%atom,method_type=method,relativistic=reltyp,error=error) 00232 00233 IF(atom_consistent_method(method,state%multiplicity)) THEN 00234 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%METHOD_INFO",extension=".log",error=error) 00235 CALL atom_print_method(atom_info(in,im)%atom,iw,error) 00236 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%METHOD_INFO",error=error) 00237 00238 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%POTENTIAL",extension=".log",error=error) 00239 IF ( pp_calc ) THEN 00240 IF(iw > 0) CALL atom_print_potential(p_pot,iw,error) 00241 ELSE 00242 IF(iw > 0) CALL atom_print_potential(ae_pot,iw,error) 00243 END IF 00244 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%POTENTIAL",error=error) 00245 END IF 00246 00247 ! calculate integrals 00248 IF ( pp_calc ) THEN 00249 ! general integrals 00250 CALL atom_int_setup(pp_int,pp_basis,& 00251 potential=p_pot,eri_coulomb=eri_c,eri_exchange=eri_e,error=error) 00252 ! potential 00253 CALL atom_ppint_setup(pp_int,pp_basis,potential=p_pot,error=error) 00254 ! 00255 NULLIFY(pp_int%tzora,pp_int%hdkh) 00256 ! 00257 CALL set_atom(atom_info(in,im)%atom,basis=pp_basis,integrals=pp_int,potential=p_pot,error=error) 00258 state%maxn_calc(:) = MIN( state%maxn_calc(:), pp_basis%nbas(:) ) 00259 CPPostcondition(ALL(state%maxn_calc(:) >= state%maxn_occ), cp_failure_level, routineP, error, failure) 00260 had_pp = .TRUE. 00261 ELSE 00262 ! general integrals 00263 CALL atom_int_setup(ae_int,ae_basis,potential=ae_pot,& 00264 eri_coulomb=eri_c,eri_exchange=eri_e,error=error) 00265 ! potential 00266 CALL atom_ppint_setup(ae_int,ae_basis,potential=ae_pot,error=error) 00267 ! relativistic correction terms 00268 CALL atom_relint_setup(ae_int,ae_basis,reltyp,zcore=REAL(zval,dp),error=error) 00269 ! 00270 CALL set_atom(atom_info(in,im)%atom,basis=ae_basis,integrals=ae_int,potential=ae_pot,error=error) 00271 state%maxn_calc(:) = MIN( state%maxn_calc(:), ae_basis%nbas(:) ) 00272 CPPostcondition(ALL(state%maxn_calc(:) >= state%maxn_occ), cp_failure_level, routineP, error, failure) 00273 had_ae = .TRUE. 00274 END IF 00275 00276 CALL set_atom(atom_info(in,im)%atom,state=state,error=error) 00277 00278 CALL set_atom(atom_info(in,im)%atom,coulomb_integral_type=do_eric,& 00279 exchange_integral_type=do_erie,error=error) 00280 00281 NULLIFY(orbitals) 00282 mo = MAXVAL(state%maxn_calc) 00283 mb = MAXVAL(atom_info(in,im)%atom%basis%nbas) 00284 CALL create_atom_orbs(orbitals,mb,mo,error) 00285 CALL set_atom(atom_info(in,im)%atom,orbitals=orbitals,error=error) 00286 00287 IF(atom_consistent_method(method,state%multiplicity)) THEN 00288 !Calculate the electronic structure 00289 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%SCF_INFO",extension=".log",error=error) 00290 CALL calculate_atom(atom_info(in,im)%atom,iw,error=error) 00291 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%SCF_INFO",error=error) 00292 00293 ! Print out the orbitals if requested 00294 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%ORBITALS",extension=".log",error=error) 00295 IF (iw > 0) THEN 00296 CALL atom_print_orbitals(atom_info(in,im)%atom,iw,error=error) 00297 END IF 00298 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%ORBITALS",error=error) 00299 00300 ! perform a fit of the total electronic density 00301 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%FIT_DENSITY",extension=".log",error=error) 00302 IF (iw>0) THEN 00303 CALL section_vals_val_get(atom_section,"PRINT%FIT_DENSITY%NUM_GTO", i_val=num_gto, error=error) 00304 powell_section => section_vals_get_subs_vals(atom_section,"POWELL",error=error) 00305 CALL atom_fit_density(atom_info(in,im)%atom,0,num_gto,iw,powell_section=powell_section,error=error) 00306 END IF 00307 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%FIT_DENSITY",error=error) 00308 00309 ! generate a response basis 00310 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%RESPONSE_BASIS",extension=".log",error=error) 00311 IF (iw>0) THEN 00312 CALL section_vals_val_get(atom_section,"PRINT%RESPONSE_BASIS%DELTA_CHARGE", r_val=delta, error=error) 00313 CALL section_vals_val_get(atom_section,"PRINT%RESPONSE_BASIS%DERIVATIVES", i_val=nder, error=error) 00314 CALL atom_response_basis(atom_info(in,im)%atom,delta,nder,iw,error) 00315 END IF 00316 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%RESPONSE_BASIS",error=error) 00317 00318 ! generate a UPF file 00319 iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%UPF_FILE",extension=".upf",& 00320 file_position="REWIND",error=error) 00321 IF (iw>0) THEN 00322 CALL atom_write_upf(atom_info(in,im)%atom,iw,error) 00323 END IF 00324 CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%UPF_FILE",error=error) 00325 END IF 00326 00327 END DO 00328 END DO 00329 00330 ! clean up 00331 IF ( had_ae ) THEN 00332 CALL atom_int_release(ae_int,error) 00333 CALL atom_ppint_release(ae_int,error) 00334 CALL atom_relint_release(ae_int,error) 00335 END IF 00336 IF ( had_pp ) THEN 00337 CALL atom_int_release(pp_int,error) 00338 CALL atom_ppint_release(pp_int,error) 00339 CALL atom_relint_release(pp_int,error) 00340 END IF 00341 CALL release_atom_basis(ae_basis,error) 00342 CALL release_atom_basis(pp_basis,error) 00343 00344 CALL release_atom_potential(p_pot,error) 00345 CALL release_atom_potential(ae_pot,error) 00346 00347 DO in = 1, n_rep 00348 DO im = 1, n_meth 00349 CALL release_atom_type(atom_info(in,im)%atom,error) 00350 END DO 00351 END DO 00352 DEALLOCATE(atom_info,STAT=ierr) 00353 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00354 00355 DEALLOCATE(ae_pot,p_pot,ae_basis,pp_basis,ae_int,pp_int,STAT=ierr) 00356 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00357 00358 CALL timestop(handle) 00359 00360 END SUBROUTINE atom_energy_opt 00361 00362 ! ***************************************************************************** 00363 SUBROUTINE atom_response_basis(atom,delta,nder,iw,error) 00364 00365 TYPE(atom_type), POINTER :: atom 00366 REAL(KIND=dp), INTENT(IN) :: delta 00367 INTEGER, INTENT(IN) :: nder, iw 00368 TYPE(cp_error_type), INTENT(INOUT) :: error 00369 00370 CHARACTER(len=*), PARAMETER :: routineN = 'atom_response_basis', 00371 routineP = moduleN//':'//routineN 00372 00373 INTEGER :: i, ider, ierr, j, k, l, 00374 lhomo, m, n, nhomo, s1, s2 00375 LOGICAL :: failure = .FALSE. 00376 REAL(KIND=dp) :: dene, emax, expzet, fhomo, o, 00377 prefac, zeta 00378 REAL(KIND=dp), ALLOCATABLE, 00379 DIMENSION(:, :) :: amat 00380 REAL(KIND=dp), ALLOCATABLE, 00381 DIMENSION(:, :, :) :: rbasis 00382 REAL(KIND=dp), ALLOCATABLE, 00383 DIMENSION(:, :, :, :) :: wfn 00384 REAL(KIND=dp), DIMENSION(:, :, :), 00385 POINTER :: ovlp 00386 TYPE(atom_state), POINTER :: state 00387 00388 WRITE(iw,'(/," ",79("*"),/,T34,A,/," ",79("*"))') "RESPONSE BASIS" 00389 00390 state => atom%state 00391 ovlp => atom%integrals%ovlp 00392 00393 ! find HOMO 00394 lhomo = -1 00395 nhomo = -1 00396 emax=-HUGE(1._dp) 00397 DO l=0,state%maxl_occ 00398 DO i=1,state%maxn_occ(l) 00399 IF (atom%orbitals%ener(i,l) > emax) THEN 00400 lhomo=l 00401 nhomo=i 00402 emax=atom%orbitals%ener(i,l) 00403 fhomo=state%occupation(l,i) 00404 END IF 00405 END DO 00406 END DO 00407 00408 s1=SIZE(atom%orbitals%wfn,1) 00409 s2=SIZE(atom%orbitals%wfn,2) 00410 ALLOCATE(wfn(s1,s2,0:3,-nder:nder),STAT=ierr) 00411 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00412 s2=MAXVAL(state%maxn_occ)+nder 00413 ALLOCATE(rbasis(s1,s2,0:3),STAT=ierr) 00414 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00415 rbasis = 0._dp 00416 00417 DO ider=-nder,nder 00418 dene = REAL(ider,KIND=dp)*delta 00419 CPPostcondition(fhomo>ABS(dene), cp_failure_level, routineP, error, failure) 00420 state%occupation(lhomo,nhomo)=fhomo+dene 00421 CALL calculate_atom(atom,iw=0,noguess=.TRUE.,error=error) 00422 wfn(:,:,:,ider) = atom%orbitals%wfn 00423 state%occupation(lhomo,nhomo)=fhomo 00424 END DO 00425 00426 DO l=0,state%maxl_occ 00427 ! occupied states 00428 DO i=1,MAX(state%maxn_occ(l),1) 00429 rbasis(:,i,l) = wfn(:,i,l,0) 00430 END DO 00431 ! differentiation 00432 DO ider=1,nder 00433 i=MAX(state%maxn_occ(l),1) 00434 SELECT CASE (ider) 00435 CASE (1) 00436 rbasis(:,i+1,l) = 0.5_dp*(wfn(:,i,l,1) - wfn(:,i,l,-1))/delta 00437 CASE (2) 00438 rbasis(:,i+2,l) = 0.25_dp*(wfn(:,i,l,2) - 2._dp*wfn(:,i,l,0) + wfn(:,i,l,-2))/delta**2 00439 CASE (3) 00440 rbasis(:,i+3,l) = 0.125_dp*(wfn(:,i,l,3) - 3._dp*wfn(:,i,l,1) & 00441 + 3._dp*wfn(:,i,l,-1) - wfn(:,i,l,-3))/delta**3 00442 CASE DEFAULT 00443 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00444 END SELECT 00445 END DO 00446 00447 ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn. 00448 n = state%maxn_occ(l)+nder 00449 m = atom%basis%nbas(l) 00450 DO i=1,n 00451 DO j=1,i-1 00452 o=DOT_PRODUCT(rbasis(1:m,j,l),RESHAPE(MATMUL(ovlp(1:m,1:m,l),rbasis(1:m,i:i,l)),(/m/))) 00453 rbasis(1:m,i,l)=rbasis(1:m,i,l)-o*rbasis(1:m,j,l) 00454 ENDDO 00455 o=DOT_PRODUCT(rbasis(1:m,i,l),RESHAPE(MATMUL(ovlp(1:m,1:m,l),rbasis(1:m,i:i,l)),(/m/))) 00456 rbasis(1:m,i,l)=rbasis(1:m,i,l)/SQRT(o) 00457 ENDDO 00458 00459 ! check 00460 ALLOCATE(amat(n,n),STAT=ierr) 00461 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00462 amat(1:n,1:n) = MATMUL(TRANSPOSE(rbasis(1:m,1:n,l)),MATMUL(ovlp(1:m,1:m,l),rbasis(1:m,1:n,l))) 00463 DO i=1,n 00464 amat(i,i)=amat(i,i) - 1._dp 00465 END DO 00466 IF (MAXVAL(ABS(amat)) > 1.e-12) THEN 00467 WRITE(iw,'(A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat)) 00468 END IF 00469 DEALLOCATE(amat,STAT=ierr) 00470 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00471 00472 ! Quickstep normalization 00473 WRITE(iw,'(/,A,T30,I3)') " Angular momentum :",l 00474 00475 WRITE(iw,'(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ",& 00476 n-nder," valence ",nder," response" 00477 expzet = 0.25_dp*REAL(2*l + 3,dp) 00478 prefac = SQRT( SQRT(pi)/2._dp**(l+2)*dfac(2*l+1) ) 00479 DO i=1,m 00480 zeta = (2._dp*atom%basis%am(i,l))**expzet 00481 WRITE(iw,'(4X,F20.10,4X,15ES20.6)') atom%basis%am(i,l),((prefac*rbasis(i,k,l)/zeta),k=1,n) 00482 END DO 00483 00484 END DO 00485 00486 DEALLOCATE(wfn,rbasis,STAT=ierr) 00487 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00488 00489 WRITE(iw,'(" ",79("*"))') 00490 00491 END SUBROUTINE atom_response_basis 00492 ! ***************************************************************************** 00493 SUBROUTINE atom_write_upf(atom,iw,error) 00494 00495 TYPE(atom_type), POINTER :: atom 00496 INTEGER, INTENT(IN) :: iw 00497 TYPE(cp_error_type), INTENT(INOUT) :: error 00498 00499 CHARACTER(len=*), PARAMETER :: routineN = 'atom_write_upf', 00500 routineP = moduleN//':'//routineN 00501 00502 CHARACTER(len=default_string_length) :: string 00503 INTEGER :: i, ibeta, ierr, j, k, l, 00504 lmax, nbeta, nr, nwfc, nwfn 00505 LOGICAL :: failure, up 00506 REAL(KIND=dp) :: pf, rl, rmax 00507 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, 00508 locpot, rp 00509 REAL(KIND=dp), ALLOCATABLE, 00510 DIMENSION(:, :) :: dij 00511 TYPE(atom_gthpot_type), POINTER :: pot 00512 00513 failure = .FALSE. 00514 IF(.NOT.atom%pp_calc) RETURN 00515 IF(atom%potential%ppot_type /= GTH_PSEUDO) RETURN 00516 pot => atom%potential%gth_pot 00517 CPPostcondition(.NOT.pot%lsdpot, cp_failure_level, routineP, error, failure) 00518 00519 WRITE(iw,'(A)') '<UPF version="2.0.1">' 00520 00521 WRITE(iw,'(T4,A)') '<PP_INFO>' 00522 WRITE(iw,'(T8,A)') 'Converted from CP2K GTH format' 00523 WRITE(iw,'(T8,A)') '<PP_INPUTFILE>' 00524 CALL atom_write_pseudo_param(pot,iw,error) 00525 WRITE(iw,'(T8,A)') '</PP_INPUTFILE>' 00526 WRITE(iw,'(T4,A)') '</PP_INFO>' 00527 WRITE(iw,'(T4,A)') '<PP_HEADER' 00528 WRITE(iw,'(T8,A)') 'generated="Generated in analytical, separable form"' 00529 WRITE(iw,'(T8,A)') 'author="Goedecker/Hartwigsen/Hutter/Teter"' 00530 WRITE(iw,'(T8,A)') 'date="Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"' 00531 WRITE(iw,'(T8,A)') 'comment="This file generated by CP2K ATOM code"' 00532 CALL compose(string,"element",cval=pot%symbol) 00533 WRITE(iw,'(T8,A)') TRIM(string) 00534 WRITE(iw,'(T8,A)') 'pseudo_type="NC"' 00535 WRITE(iw,'(T8,A)') 'relativistic="no"' 00536 WRITE(iw,'(T8,A)') 'is_ultrasoft="F"' 00537 WRITE(iw,'(T8,A)') 'is_paw="F"' 00538 WRITE(iw,'(T8,A)') 'is_coulomb="F"' 00539 WRITE(iw,'(T8,A)') 'has_so="F"' 00540 WRITE(iw,'(T8,A)') 'has_wfc="F"' 00541 WRITE(iw,'(T8,A)') 'has_gipaw="F"' 00542 WRITE(iw,'(T8,A)') 'paw_as_gipaw="F"' 00543 IF(pot%nlcc) THEN 00544 WRITE(iw,'(T8,A)') 'core_correction="T"' 00545 ELSE 00546 WRITE(iw,'(T8,A)') 'core_correction="F"' 00547 END IF 00548 WRITE(iw,'(T8,A)') 'functional="DFT"' 00549 CALL compose(string,"z_valence",rval=pot%zion) 00550 WRITE(iw,'(T8,A)') TRIM(string) 00551 CALL compose(string,"total_psenergy",rval=2._dp*atom%energy%etot) 00552 WRITE(iw,'(T8,A)') TRIM(string) 00553 WRITE(iw,'(T8,A)') 'wfc_cutoff="0.0E+00"' 00554 WRITE(iw,'(T8,A)') 'rho_cutoff="0.0E+00"' 00555 lmax = -1 00556 DO l=0,3 00557 IF(pot%nl(l) > 0) lmax = l 00558 END DO 00559 CALL compose(string,"l_max",ival=lmax) 00560 WRITE(iw,'(T8,A)') TRIM(string) 00561 WRITE(iw,'(T8,A)') 'l_max_rho="0"' 00562 WRITE(iw,'(T8,A)') 'l_local="-3"' 00563 nr = atom%basis%grid%nr 00564 CALL compose(string,"mesh_size",ival=nr) 00565 WRITE(iw,'(T8,A)') TRIM(string) 00566 nwfc = SUM(atom%state%maxn_occ) 00567 CALL compose(string,"number_of_wfc",ival=nwfc) 00568 WRITE(iw,'(T8,A)') TRIM(string) 00569 nbeta = SUM(pot%nl) 00570 CALL compose(string,"number_of_proj",ival=nbeta) 00571 WRITE(iw,'(T8,A)') TRIM(string)//'/>' 00572 00573 ! Mesh 00574 up = atom%basis%grid%rad(1) < atom%basis%grid%rad(nr) 00575 WRITE(iw,'(T4,A)') '<PP_MESH' 00576 WRITE(iw,'(T8,A)') 'dx="1.E+00"' 00577 CALL compose(string,"mesh",ival=nr) 00578 WRITE(iw,'(T8,A)') TRIM(string) 00579 WRITE(iw,'(T8,A)') 'xmin="1.E+00"' 00580 rmax = MAXVAL(atom%basis%grid%rad) 00581 CALL compose(string,"rmax",rval=rmax) 00582 WRITE(iw,'(T8,A)') TRIM(string) 00583 WRITE(iw,'(T8,A)') 'zmesh="1.E+00">' 00584 WRITE(iw,'(T8,A)') '<PP_R type="real"' 00585 CALL compose(string,"size",ival=nr) 00586 WRITE(iw,'(T12,A)') TRIM(string) 00587 WRITE(iw,'(T12,A)') 'columns="4">' 00588 IF(up) THEN 00589 WRITE(iw,'(T12,4ES25.12E3)') (atom%basis%grid%rad(i),i=1,nr) 00590 ELSE 00591 WRITE(iw,'(T12,4ES25.12E3)') (atom%basis%grid%rad(i),i=nr,1,-1) 00592 END IF 00593 WRITE(iw,'(T8,A)') '</PP_R>' 00594 WRITE(iw,'(T8,A)') '<PP_RAB type="real"' 00595 CALL compose(string,"size",ival=nr) 00596 WRITE(iw,'(T12,A)') TRIM(string) 00597 WRITE(iw,'(T12,A)') 'columns="4">' 00598 IF(up) THEN 00599 WRITE(iw,'(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i),i=1,nr) 00600 ELSE 00601 WRITE(iw,'(T12,4ES25.12E3)') (atom%basis%grid%wr(i)/atom%basis%grid%rad2(i),i=nr,1,-1) 00602 END IF 00603 WRITE(iw,'(T8,A)') '</PP_RAB>' 00604 WRITE(iw,'(T4,A)') '</PP_MESH>' 00605 00606 ! NLCC 00607 IF(pot%nlcc) THEN 00608 WRITE(iw,'(T4,A)') '<PP_NLCC type="real"' 00609 CALL compose(string,"size",ival=nr) 00610 WRITE(iw,'(T8,A)') TRIM(string) 00611 WRITE(iw,'(T8,A)') 'columns="4">' 00612 ALLOCATE(corden(nr),STAT=ierr) 00613 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00614 CALL atom_core_density(corden,pot,"RHO",atom%basis%grid%rad,error) 00615 IF(up) THEN 00616 WRITE(iw,'(T8,4ES25.12E3)') (corden(i),i=1,nr) 00617 ELSE 00618 WRITE(iw,'(T8,4ES25.12E3)') (corden(i),i=nr,1,-1) 00619 END IF 00620 DEALLOCATE(corden,STAT=ierr) 00621 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00622 WRITE(iw,'(T4,A)') '</PP_NLCC>' 00623 END IF 00624 00625 ! local PP 00626 WRITE(iw,'(T4,A)') '<PP_LOCAL type="real"' 00627 CALL compose(string,"size",ival=nr) 00628 WRITE(iw,'(T8,A)') TRIM(string) 00629 WRITE(iw,'(T8,A)') 'columns="4">' 00630 ALLOCATE(locpot(nr),STAT=ierr) 00631 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00632 CALL atom_local_potential(locpot,pot,atom%basis%grid%rad,error) 00633 IF(up) THEN 00634 WRITE(iw,'(T8,4ES25.12E3)') (2.0_dp*locpot(i),i=1,nr) 00635 ELSE 00636 WRITE(iw,'(T8,4ES25.12E3)') (2.0_dp*locpot(i),i=nr,1,-1) 00637 END IF 00638 DEALLOCATE(locpot,STAT=ierr) 00639 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00640 WRITE(iw,'(T4,A)') '</PP_LOCAL>' 00641 00642 ! nonlocal PP 00643 WRITE(iw,'(T4,A)') '<PP_NONLOCAL>' 00644 ALLOCATE(rp(nr),ef(nr),beta(nr),STAT=ierr) 00645 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00646 ibeta = 0 00647 DO l=0,3 00648 IF(pot%nl(l) == 0) CYCLE 00649 rl = pot%rcnl(l) 00650 rp = atom%basis%grid%rad 00651 ef = EXP(-0.5_dp*rp*rp/(rl*rl)) 00652 DO i=1,pot%nl(l) 00653 pf = rl**(l+0.5_dp*(4._dp*i-1._dp)) 00654 j = l + 2*i - 1 00655 pf = SQRT(2._dp) / (pf * SQRT(gamma1(j))) 00656 beta = pf * rp**(l+2*i-2) * ef 00657 ibeta=ibeta+1 00658 CALL compose(string,"<PP_BETA",counter=ibeta) 00659 WRITE(iw,'(T8,A)') TRIM(string) 00660 CALL compose(string,"angular_momentum",ival=l) 00661 WRITE(iw,'(T12,A)') TRIM(string) 00662 WRITE(iw,'(T12,A)') 'type="real"' 00663 CALL compose(string,"size",ival=nr) 00664 WRITE(iw,'(T12,A)') TRIM(string) 00665 WRITE(iw,'(T12,A)') 'columns="4">' 00666 beta = beta*rp 00667 IF(up) THEN 00668 WRITE(iw,'(T12,4ES25.12E3)') (2._dp*beta(j),j=1,nr) 00669 ELSE 00670 WRITE(iw,'(T12,4ES25.12E3)') (2._dp*beta(j),j=nr,1,-1) 00671 END IF 00672 CALL compose(string,"</PP_BETA",counter=ibeta,isfinal=.TRUE.) 00673 WRITE(iw,'(T8,A)') TRIM(string) 00674 END DO 00675 END DO 00676 DEALLOCATE(rp,ef,beta,STAT=ierr) 00677 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00678 ! nonlocal PP matrix elements 00679 ALLOCATE(dij(nbeta,nbeta),STAT=ierr) 00680 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00681 dij = 0._dp 00682 DO l=0,3 00683 IF(pot%nl(l) == 0) CYCLE 00684 ibeta = SUM(pot%nl(0:l-1)) + 1 00685 i = ibeta + pot%nl(l) - 1 00686 dij(ibeta:i,ibeta:i) = pot%hnl(1:pot%nl(l),1:pot%nl(l),l) 00687 END DO 00688 WRITE(iw,'(T8,A)') '<PP_DIJ type="real"' 00689 WRITE(iw,'(T12,A)') 'columns="4">' 00690 WRITE(iw,'(T12,4ES25.12E3)') ((0.5_dp*dij(i,j),j=1,nbeta),i=1,nbeta) 00691 WRITE(iw,'(T8,A)') '</PP_DIJ>' 00692 DEALLOCATE(dij,STAT=ierr) 00693 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00694 WRITE(iw,'(T4,A)') '</PP_NONLOCAL>' 00695 00696 ! atomic wavefunctions 00697 ALLOCATE(beta(nr),STAT=ierr) 00698 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00699 WRITE(iw,'(T4,A)') '<PP_PSWFC>' 00700 nwfn=0 00701 DO l=0,3 00702 DO i=1,10 00703 IF(ABS(atom%state%occupation(l,i)) == 0._dp) CYCLE 00704 nwfn = nwfn+1 00705 CALL compose(string,"<PP_CHI",counter=nwfn) 00706 WRITE(iw,'(T8,A)') TRIM(string) 00707 CALL compose(string,"l",ival=l) 00708 WRITE(iw,'(T12,A)') TRIM(string) 00709 CALL compose(string,"occupation",rval=atom%state%occupation(l,i)) 00710 WRITE(iw,'(T12,A)') TRIM(string) 00711 CALL compose(string,"pseudo_energy",rval=2._dp*atom%orbitals%ener(i,l)) 00712 WRITE(iw,'(T12,A)') TRIM(string) 00713 WRITE(iw,'(T12,A)') 'columns="4">' 00714 beta = 0._dp 00715 DO k=1,atom%basis%nbas(l) 00716 beta(:) = beta(:) + atom%orbitals%wfn(k,i,l)*atom%basis%bf(:,k,l) 00717 END DO 00718 beta = beta*atom%basis%grid%rad 00719 IF(up) THEN 00720 WRITE(iw,'(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j),j=1,nr) 00721 ELSE 00722 WRITE(iw,'(T12,4ES25.12E3)') (beta(j)*atom%basis%grid%rad(j),j=nr,1,-1) 00723 END IF 00724 CALL compose(string,"</PP_CHI",counter=nwfn,isfinal=.TRUE.) 00725 WRITE(iw,'(T8,A)') TRIM(string) 00726 END DO 00727 END DO 00728 WRITE(iw,'(T4,A)') '</PP_PSWFC>' 00729 DEALLOCATE(beta,STAT=ierr) 00730 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00731 00732 ! atomic charge 00733 ALLOCATE(dens(nr),STAT=ierr) 00734 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00735 WRITE(iw,'(T4,A)') '<PP_RHOATOM type="real"' 00736 CALL compose(string,"size",ival=nr) 00737 WRITE(iw,'(T8,A)') TRIM(string) 00738 WRITE(iw,'(T8,A)') 'columns="4">' 00739 CALL atom_density(dens,atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,& 00740 "RHO",atom%basis%grid%rad,error) 00741 IF(up) THEN 00742 WRITE(iw,'(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j),j=1,nr) 00743 ELSE 00744 WRITE(iw,'(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j),j=nr,1,-1) 00745 END IF 00746 WRITE(iw,'(T4,A)') '</PP_RHOATOM>' 00747 DEALLOCATE(dens,STAT=ierr) 00748 CPPostcondition(ierr==0,cp_failure_level,routineP,error,failure) 00749 00750 WRITE(iw,'(A)') '</UPF>' 00751 00752 END SUBROUTINE atom_write_upf 00753 SUBROUTINE compose(string,tag,counter,rval,ival,cval,isfinal) 00754 CHARACTER(len=*) :: string, tag 00755 INTEGER, OPTIONAL :: counter 00756 REAL(KIND=dp), OPTIONAL :: rval 00757 INTEGER, OPTIONAL :: ival 00758 CHARACTER(len=*), OPTIONAL :: cval 00759 LOGICAL, OPTIONAL :: isfinal 00760 00761 CHARACTER(len=default_string_length) :: str 00762 LOGICAL :: fin 00763 00764 IF(PRESENT(counter)) THEN 00765 WRITE(str,"(I12)") counter 00766 ELSEIF(PRESENT(rval)) THEN 00767 WRITE(str,"(G18.8)") rval 00768 ELSEIF(PRESENT(ival)) THEN 00769 WRITE(str,"(I12)") ival 00770 ELSEIF(PRESENT(cval)) THEN 00771 WRITE(str,"(A)") TRIM(ADJUSTL(cval)) 00772 ELSE 00773 WRITE(str,"(A)") "" 00774 END IF 00775 fin = .FALSE. 00776 IF(PRESENT(isfinal)) fin=isfinal 00777 IF(PRESENT(counter)) THEN 00778 IF(fin) THEN 00779 WRITE(string,"(A,A1,A,A1)") TRIM(ADJUSTL(tag)),'.',TRIM(ADJUSTL(str)),'>' 00780 ELSE 00781 WRITE(string,"(A,A1,A)") TRIM(ADJUSTL(tag)),'.',TRIM(ADJUSTL(str)) 00782 END IF 00783 ELSE 00784 IF(fin) THEN 00785 WRITE(string,"(A,A2,A,A2)") TRIM(ADJUSTL(tag)),'="',TRIM(ADJUSTL(str)),'>"' 00786 ELSE 00787 WRITE(string,"(A,A2,A,A1)") TRIM(ADJUSTL(tag)),'="',TRIM(ADJUSTL(str)),'"' 00788 END IF 00789 END IF 00790 00791 END SUBROUTINE compose 00792 ! ***************************************************************************** 00793 00794 END MODULE atom_energy
1.7.3