CP2K 2.4 (Revision 12889)

atom_types.f90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------------!
00002 !   CP2K: A general program to perform molecular dynamics simulations         !
00003 !   Copyright (C) 2000 - 2013  CP2K developers group                          !
00004 !-----------------------------------------------------------------------------!
00005 
00006 ! *****************************************************************************
00013 MODULE atom_types
00014   USE cp_linked_list_val,              ONLY: cp_sll_val_next,&
00015                                              cp_sll_val_type
00016   USE cp_parser_methods,               ONLY: parser_get_next_line,&
00017                                              parser_get_object,&
00018                                              parser_search_string,&
00019                                              parser_test_next_token
00020   USE cp_parser_types,                 ONLY: cp_parser_type,&
00021                                              parser_create,&
00022                                              parser_release
00023   USE f77_blas
00024   USE input_constants,                 ONLY: contracted_gto,&
00025                                              gaussian,&
00026                                              geometrical_gto,&
00027                                              numerical,&
00028                                              slater
00029   USE input_section_types,             ONLY: section_vals_get,&
00030                                              section_vals_get_subs_vals,&
00031                                              section_vals_list_get,&
00032                                              section_vals_type,&
00033                                              section_vals_val_get
00034   USE input_val_types,                 ONLY: val_get,&
00035                                              val_type
00036   USE kinds,                           ONLY: default_string_length,&
00037                                              dp
00038   USE mathconstants,                   ONLY: dfac,&
00039                                              fac,&
00040                                              pi
00041   USE periodic_table,                  ONLY: ptable
00042   USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
00043                                              create_grid_atom,&
00044                                              deallocate_grid_atom,&
00045                                              grid_atom_type
00046   USE string_utilities,                ONLY: remove_word,&
00047                                              uppercase
00048 #include "cp_common_uses.h"
00049 
00050   IMPLICIT NONE
00051 
00052   PRIVATE
00053 
00054   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_types'
00055 
00056   INTEGER, PARAMETER                                 :: GTO_BASIS=100,
00057                                                         CGTO_BASIS=101,
00058                                                         STO_BASIS=102,
00059                                                         NUM_BASIS=103
00060   INTEGER, PARAMETER                                 :: NO_PSEUDO=0,
00061                                                         GTH_PSEUDO=1
00062 
00065   TYPE atom_basis_type
00066     INTEGER                                       :: basis_type
00067     INTEGER, DIMENSION(0:3)                       :: nbas
00068     INTEGER, DIMENSION(0:3)                       :: nprim
00069     REAL(KIND=dp),DIMENSION(:,:),POINTER          :: am         !GTO exponents
00070     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: cm         !Contraction coeffs
00071     REAL(KIND=dp),DIMENSION(:,:),POINTER          :: as         !STO exponents
00072     INTEGER,DIMENSION(:,:),POINTER                :: ns         !STO n-quantum numbers
00073     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: bf         !num. bsf
00074     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: dbf        !derivatives (num)
00075     REAL(KIND=dp)                                 :: eps_eig
00076     TYPE(grid_atom_type),POINTER                  :: grid
00077     LOGICAL                                       :: geometrical
00078     REAL(KIND=dp)                                 :: aval, cval
00079     INTEGER, DIMENSION(0:3)                       :: start
00080   END TYPE atom_basis_type
00081 
00084   TYPE atom_gthpot_type
00085     CHARACTER (LEN=2)                             :: symbol
00086     CHARACTER (LEN=80)                            :: pname
00087     INTEGER, DIMENSION(0:3)                       :: econf
00088     REAL(dp)                                      :: zion
00089     REAL(dp)                                      :: rc
00090     INTEGER                                       :: ncl
00091     REAL(dp), DIMENSION(5)                        :: cl
00092     INTEGER, DIMENSION(0:3)                       :: nl
00093     REAL(dp), DIMENSION(0:3)                      :: rcnl
00094     REAL(dp), DIMENSION(4,4,0:3)                  :: hnl
00095     ! type extensions
00096     ! NLCC
00097     LOGICAL                                    :: nlcc
00098     INTEGER                                    :: nexp_nlcc
00099     REAL(KIND = dp), DIMENSION(10)             :: alpha_nlcc
00100     INTEGER, DIMENSION(10)                     :: nct_nlcc
00101     REAL(KIND = dp), DIMENSION(4,10)           :: cval_nlcc
00102     ! LSD potential
00103     LOGICAL                                    :: lsdpot
00104     INTEGER                                    :: nexp_lsd
00105     REAL(KIND = dp), DIMENSION(10)             :: alpha_lsd
00106     INTEGER, DIMENSION(10)                     :: nct_lsd
00107     REAL(KIND = dp), DIMENSION(4,10)           :: cval_lsd
00108     ! extended local potential
00109     LOGICAL                                    :: lpotextended
00110     INTEGER                                    :: nexp_lpot
00111     REAL(KIND = dp), DIMENSION(10)             :: alpha_lpot
00112     INTEGER, DIMENSION(10)                     :: nct_lpot
00113     REAL(KIND = dp), DIMENSION(4,10)           :: cval_lpot
00114   END TYPE atom_gthpot_type
00115 
00116   TYPE atom_potential_type
00117     INTEGER                                       :: ppot_type
00118     LOGICAL                                       :: confinement
00119     REAL(dp)                                      :: acon
00120     REAL(dp)                                      :: rcon
00121     INTEGER                                       :: ncon
00122     TYPE(atom_gthpot_type)                        :: gth_pot
00123   END TYPE atom_potential_type
00124 
00127   TYPE atom_state
00128     REAL(KIND=dp),DIMENSION(0:3,10)               :: occ
00129     REAL(KIND=dp),DIMENSION(0:3,10)               :: core
00130     REAL(KIND=dp),DIMENSION(0:3,10)               :: occupation
00131     INTEGER                                       :: maxl_occ
00132     INTEGER,DIMENSION(0:3)                        :: maxn_occ
00133     INTEGER                                       :: maxl_calc
00134     INTEGER,DIMENSION(0:3)                        :: maxn_calc
00135     INTEGER                                       :: multiplicity
00136     REAL(KIND=dp),DIMENSION(0:3,10)               :: occa, occb
00137   END TYPE atom_state
00138 
00141   TYPE eri
00142     REAL(KIND=dp),DIMENSION(:,:),POINTER          :: int
00143   END TYPE eri
00144 
00145   TYPE atom_integrals
00146     INTEGER                                       :: status=0
00147     INTEGER                                       :: ppstat=0
00148     LOGICAL                                       :: eri_coulomb
00149     LOGICAL                                       :: eri_exchange
00150     LOGICAL                                       :: all_nu
00151     INTEGER, DIMENSION(0:3)                       :: n, nne
00152     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: ovlp, kin, core, clsd
00153     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: utrans, uptrans
00154     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: hnl
00155     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: conf
00156     TYPE(eri),DIMENSION(20)                       :: ceri
00157     TYPE(eri),DIMENSION(20)                       :: eeri
00158     INTEGER                                       :: dkhstat=0
00159     INTEGER                                       :: zorastat=0
00160     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: tzora
00161     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: hdkh
00162   END TYPE atom_integrals
00163 
00166   TYPE atom_orbitals
00167     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: wfn,wfna,wfnb
00168     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: pmat,pmata,pmatb
00169     REAL(KIND=dp),DIMENSION(:,:),POINTER          :: ener,enera,enerb
00170     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: refene,refchg,refnod
00171     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: wrefene,wrefchg,wrefnod
00172     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: crefene,crefchg,crefnod
00173     REAL(KIND=dp),DIMENSION(:,:),POINTER          :: wpsir0
00174     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: rcmax
00175   END TYPE atom_orbitals
00176 
00179   TYPE opmat_type
00180     INTEGER, DIMENSION(0:3)                       :: n
00181     REAL(KIND=dp),DIMENSION(:,:,:),POINTER        :: op
00182   END TYPE opmat_type
00183 
00186   TYPE opgrid_type
00187     REAL(KIND=dp),DIMENSION(:),POINTER            :: op
00188     TYPE(grid_atom_type),POINTER                  :: grid
00189   END TYPE opgrid_type
00190 
00193   TYPE atom_energy_type
00194     REAL(KIND=dp)                                 :: etot
00195     REAL(KIND=dp)                                 :: eband
00196     REAL(KIND=dp)                                 :: ekin
00197     REAL(KIND=dp)                                 :: epot
00198     REAL(KIND=dp)                                 :: ecore
00199     REAL(KIND=dp)                                 :: elsd
00200     REAL(KIND=dp)                                 :: epseudo
00201     REAL(KIND=dp)                                 :: eploc
00202     REAL(KIND=dp)                                 :: epnl
00203     REAL(KIND=dp)                                 :: exc
00204     REAL(KIND=dp)                                 :: ecoulomb
00205     REAL(KIND=dp)                                 :: eexchange
00206     REAL(KIND=dp)                                 :: econfinement
00207   END TYPE atom_energy_type
00208 
00211   TYPE atom_optimization_type
00212     REAL(KIND=dp)                                 :: damping
00213     REAL(KIND=dp)                                 :: eps_scf
00214     REAL(KIND=dp)                                 :: eps_diis
00215     INTEGER                                       :: max_iter
00216     INTEGER                                       :: n_diis
00217   END TYPE atom_optimization_type
00218 
00221   TYPE atom_type
00222     INTEGER                                       :: z
00223     INTEGER                                       :: zcore
00224     LOGICAL                                       :: pp_calc
00225     INTEGER                                       :: method_type
00226     INTEGER                                       :: relativistic
00227     INTEGER                                       :: coulomb_integral_type
00228     INTEGER                                       :: exchange_integral_type
00229     REAL(KIND=dp)                                 :: weight
00230     TYPE(atom_basis_type), POINTER                :: basis
00231     TYPE(atom_potential_type), POINTER            :: potential
00232     TYPE(atom_state), POINTER                     :: state
00233     TYPE(atom_integrals), POINTER                 :: integrals
00234     TYPE(atom_orbitals), POINTER                  :: orbitals
00235     TYPE(atom_energy_type)                        :: energy
00236     TYPE(atom_optimization_type)                  :: optimization
00237     TYPE(section_vals_type), POINTER              :: xc_section
00238   END TYPE atom_type
00239 ! *****************************************************************************
00240   TYPE atom_p_type
00241     TYPE(atom_type), POINTER                      :: atom
00242   END TYPE atom_p_type
00243 
00244   PUBLIC :: atom_p_type, atom_type, atom_basis_type, atom_state, atom_integrals
00245   PUBLIC :: atom_orbitals, atom_energy_type, eri, atom_potential_type, atom_gthpot_type
00246   PUBLIC :: atom_optimization_type
00247   PUBLIC :: create_atom_type, release_atom_type, set_atom
00248   PUBLIC :: create_atom_orbs, release_atom_orbs
00249   PUBLIC :: init_atom_basis, release_atom_basis
00250   PUBLIC :: init_atom_potential, release_atom_potential
00251   PUBLIC :: read_atom_opt_section
00252   PUBLIC :: Clementi_geobas
00253   PUBLIC :: GTO_BASIS, CGTO_BASIS, STO_BASIS, NUM_BASIS
00254   PUBLIC :: opmat_type, create_opmat, release_opmat
00255   PUBLIC :: opgrid_type, create_opgrid, release_opgrid
00256   PUBLIC :: NO_PSEUDO, GTH_PSEUDO
00257 
00258 ! *****************************************************************************
00259 
00260 CONTAINS
00261 
00262 ! *****************************************************************************
00286   SUBROUTINE init_atom_basis(basis,basis_section,zval,btyp,error)
00287     TYPE(atom_basis_type), INTENT(INOUT)     :: basis
00288     TYPE(section_vals_type), POINTER         :: basis_section
00289     INTEGER, INTENT(IN)                      :: zval
00290     CHARACTER(LEN=2)                         :: btyp
00291     TYPE(cp_error_type), INTENT(inout)       :: error
00292 
00293     CHARACTER(len=*), PARAMETER :: routineN = 'init_atom_basis', 
00294       routineP = moduleN//':'//routineN
00295     INTEGER, PARAMETER                       :: nua = 40, nup = 16
00296     REAL(KIND=dp), DIMENSION(nua), PARAMETER :: ugbs = (/0.007299_dp,
00297       0.013705_dp, 0.025733_dp, 0.048316_dp, 0.090718_dp, 0.170333_dp,
00298       0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp, 3.974902_dp,
00299       7.463317_dp,14.013204_dp, 26.311339_dp, 49.402449_dp, 92.758561_dp,
00300       174.164456_dp, 327.013024_dp,614.003114_dp, 1152.858743_dp,
00301       2164.619772_dp, 4064.312984_dp, 7631.197056_dp,14328.416324_dp,
00302       26903.186074_dp, 50513.706789_dp, 94845.070265_dp, 178082.107320_dp,
00303       334368.848683_dp, 627814.487663_dp, 1178791.123851_dp, 2213310.684886_dp
00304       ,4155735.557141_dp,7802853.046713_dp, 14650719.428954_dp,
00305       27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp,
00306       182087882.613702_dp, 341890134.751331_dp /)
00307 
00308     CHARACTER(LEN=default_string_length)     :: basis_fn, basis_name
00309     INTEGER                                  :: basistype, i, ierr, j, k, l, 
00310                                                 ll, m, ngp, nl, nr, nu, 
00311                                                 quadtype
00312     INTEGER, DIMENSION(0:3)                  :: starti
00313     INTEGER, DIMENSION(:), POINTER           :: nqm, num_gto, num_slater, 
00314                                                 sindex
00315     LOGICAL                                  :: failure
00316     REAL(KIND=dp)                            :: al, amax, aval, cval, ear, 
00317                                                 pf, rk
00318     REAL(KIND=dp), DIMENSION(:), POINTER     :: expo
00319     TYPE(section_vals_type), POINTER         :: gto_basis_section
00320 
00321 !   btyp = AE : standard all-electron basis
00322 !   btyp = PP : standard pseudopotential basis
00323 !   btyp = AA : high accuracy all-electron basis
00324 !   btyp = AP : high accuracy pseudopotential basis
00325 
00326     failure = .FALSE.
00327     NULLIFY(basis%am,basis%cm,basis%as,basis%ns,basis%bf,basis%dbf)
00328     ! get information on quadrature type and number of grid points
00329     ! allocate and initialize the atomic grid
00330     CALL allocate_grid_atom(basis%grid,error)
00331     CALL section_vals_val_get(basis_section,"QUADRATURE",i_val=quadtype,error=error)
00332     CALL section_vals_val_get(basis_section,"GRID_POINTS",i_val=ngp,error=error)
00333     CALL cp_assert(ngp > 0,cp_failure_level,cp_assertion_failed,routineP,&
00334                    "# point radial grid < 0",error,failure)
00335     CALL create_grid_atom(basis%grid,ngp,1,1,0,quadtype)
00336     basis%grid%nr = ngp
00337     basis%geometrical = .FALSE.
00338     basis%aval  = 0._dp
00339     basis%cval  = 0._dp
00340     basis%start = 0
00341 
00342     CALL section_vals_val_get(basis_section,"BASIS_TYPE",i_val=basistype,error=error)
00343     CALL section_vals_val_get(basis_section,"EPS_EIGENVALUE",r_val=basis%eps_eig,error=error)
00344     SELECT CASE (basistype)
00345       CASE DEFAULT
00346          CPAssert(.FALSE.,cp_failure_level,routineP,error,failure)
00347       CASE (gaussian)
00348          basis%basis_type = GTO_BASIS
00349          NULLIFY(num_gto)
00350          CALL section_vals_val_get(basis_section,"NUM_GTO",i_vals=num_gto,error=error)
00351          IF ( num_gto(1) < 1 ) THEN
00352            ! use default basis
00353            IF ( btyp == "AE" ) THEN
00354              nu = nua
00355            ELSEIF ( btyp == "PP" ) THEN
00356              nu = nup
00357            ELSE
00358              nu = nua
00359            ENDIF
00360            basis%nbas = nu
00361            basis%nprim = nu
00362            ALLOCATE (basis%am(nu,0:3),STAT=ierr)
00363            CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00364            basis%am(1:nu,0) = ugbs(1:nu)
00365            basis%am(1:nu,1) = ugbs(1:nu)
00366            basis%am(1:nu,2) = ugbs(1:nu)
00367            basis%am(1:nu,3) = ugbs(1:nu)
00368          ELSE
00369            basis%nbas = 0
00370            DO i=1,SIZE(num_gto)
00371              basis%nbas(i-1) = num_gto(i)
00372            END DO
00373            basis%nprim = basis%nbas
00374            m  = MAXVAL(basis%nbas)
00375            ALLOCATE (basis%am(m,0:3),STAT=ierr)
00376            CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00377            basis%am = 0._dp
00378            IF ( basis%nbas(0) > 0 ) THEN
00379              NULLIFY(expo)
00380              CALL section_vals_val_get(basis_section,"S_EXPONENTS",r_vals=expo,error=error)
00381              CPPostcondition(SIZE(expo)>=basis%nbas(0), cp_failure_level, routineP, error, failure)
00382              DO i=1,basis%nbas(0)
00383                basis%am(i,0) = expo(i)
00384              END DO
00385            END IF
00386            IF ( basis%nbas(1) > 0 ) THEN
00387              NULLIFY(expo)
00388              CALL section_vals_val_get(basis_section,"P_EXPONENTS",r_vals=expo,error=error)
00389              CPPostcondition(SIZE(expo)>=basis%nbas(1), cp_failure_level, routineP, error, failure)
00390              DO i=1,basis%nbas(1)
00391                basis%am(i,1) = expo(i)
00392              END DO
00393            END IF
00394            IF ( basis%nbas(2) > 0 ) THEN
00395              NULLIFY(expo)
00396              CALL section_vals_val_get(basis_section,"D_EXPONENTS",r_vals=expo,error=error)
00397              CPPostcondition(SIZE(expo)>=basis%nbas(2), cp_failure_level, routineP, error, failure)
00398              DO i=1,basis%nbas(2)
00399                basis%am(i,2) = expo(i)
00400              END DO
00401            END IF
00402            IF ( basis%nbas(3) > 0 ) THEN
00403              NULLIFY(expo)
00404              CALL section_vals_val_get(basis_section,"F_EXPONENTS",r_vals=expo,error=error)
00405              CPPostcondition(SIZE(expo)>=basis%nbas(3), cp_failure_level, routineP, error, failure)
00406              DO i=1,basis%nbas(3)
00407                basis%am(i,3) = expo(i)
00408              END DO
00409            END IF
00410          END IF
00411          ! initialize basis function on a radial grid
00412          nr = basis%grid%nr
00413          m  = MAXVAL(basis%nbas)
00414          ALLOCATE (basis%bf(nr,m,0:3),STAT=ierr)
00415          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00416          ALLOCATE (basis%dbf(nr,m,0:3),STAT=ierr)
00417          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00418          basis%bf =  0._dp
00419          basis%dbf = 0._dp
00420          DO l=0,3
00421            DO i=1,basis%nbas(l)
00422              al  = basis%am(i,l)
00423              DO k=1,nr
00424                rk  = basis%grid%rad(k)
00425                ear = EXP(-al*basis%grid%rad(k)**2)
00426                basis%bf(k,i,l) = rk**l * ear
00427                basis%dbf(k,i,l) = ( REAL(l,dp)*rk**(l-1) - 2._dp*al*rk**(l+1) ) * ear
00428              END DO
00429            END DO
00430          END DO
00431       CASE (geometrical_gto)
00432          basis%basis_type = GTO_BASIS
00433          NULLIFY(num_gto)
00434          CALL section_vals_val_get(basis_section,"NUM_GTO",i_vals=num_gto,error=error)
00435          IF ( num_gto(1) < 1 ) THEN
00436            IF ( btyp == "AE" ) THEN
00437              ! use the Clementi extra large basis
00438              CALL Clementi_geobas(zval,cval,aval,basis%nbas,starti,error)
00439            ELSEIF ( btyp == "PP" ) THEN
00440              ! use the Clementi extra large basis
00441              CALL Clementi_geobas(zval,cval,aval,basis%nbas,starti,error)
00442            ELSEIF ( btyp == "AA" ) THEN
00443              CALL Clementi_geobas(zval,cval,aval,basis%nbas,starti,error)
00444              amax = cval**(basis%nbas(0)-1)
00445              basis%nbas(0) = NINT((LOG(amax)/LOG(1.6_dp)))
00446              cval = 1.6_dp
00447              starti = 0
00448              basis%nbas(1) = basis%nbas(0) - 4
00449              basis%nbas(2) = basis%nbas(0) - 8
00450              basis%nbas(3) = basis%nbas(0) -12
00451            ELSEIF ( btyp == "AP" ) THEN
00452              CALL Clementi_geobas(zval,cval,aval,basis%nbas,starti,error)
00453              amax = 500._dp/aval
00454              basis%nbas = NINT((LOG(amax)/LOG(1.6_dp)))
00455              cval = 1.6_dp
00456              starti = 0
00457            ELSE
00458              ! use the Clementi extra large basis
00459              CALL Clementi_geobas(zval,cval,aval,basis%nbas,starti,error)
00460            ENDIF
00461            basis%nprim = basis%nbas
00462          ELSE
00463            basis%nbas = 0
00464            DO i=1,SIZE(num_gto)
00465              basis%nbas(i-1) = num_gto(i)
00466            END DO
00467            basis%nprim = basis%nbas
00468            NULLIFY(sindex)
00469            CALL section_vals_val_get(basis_section,"START_INDEX",i_vals=sindex,error=error)
00470            starti = 0
00471            DO i=1,SIZE(sindex)
00472              starti(i-1) = sindex(i)
00473              CPPostcondition(sindex(i)>=0, cp_failure_level, routineP, error, failure)
00474            END DO
00475            CALL section_vals_val_get(basis_section,"GEOMETRICAL_FACTOR",r_val=cval,error=error)
00476            CALL section_vals_val_get(basis_section,"GEO_START_VALUE",r_val=aval,error=error)
00477          END IF
00478          m  = MAXVAL(basis%nbas)
00479          ALLOCATE (basis%am(m,0:3),STAT=ierr)
00480          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00481          basis%am = 0._dp
00482          DO l=0,3
00483            DO i=1,basis%nbas(l)
00484              ll = i - 1 + starti(l)
00485              basis%am(i,l) = aval * cval**(ll)
00486            END DO
00487          END DO
00488 
00489          basis%geometrical = .TRUE.
00490          basis%aval = aval
00491          basis%cval = cval
00492          basis%start = starti
00493 
00494          ! initialize basis function on a radial grid
00495          nr = basis%grid%nr
00496          m  = MAXVAL(basis%nbas)
00497          ALLOCATE (basis%bf(nr,m,0:3),STAT=ierr)
00498          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00499          ALLOCATE (basis%dbf(nr,m,0:3),STAT=ierr)
00500          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00501          basis%bf =  0._dp
00502          basis%dbf = 0._dp
00503          DO l=0,3
00504            DO i=1,basis%nbas(l)
00505              al  = basis%am(i,l)
00506              DO k=1,nr
00507                rk  = basis%grid%rad(k)
00508                ear = EXP(-al*basis%grid%rad(k)**2)
00509                basis%bf(k,i,l) = rk**l * ear
00510                basis%dbf(k,i,l) = ( REAL(l,dp)*rk**(l-1) - 2._dp*al*rk**(l+1) ) * ear
00511              END DO
00512            END DO
00513          END DO
00514       CASE (contracted_gto)
00515          basis%basis_type = CGTO_BASIS
00516          CALL section_vals_val_get(basis_section,"BASIS_SET_FILE_NAME",c_val=basis_fn,error=error)
00517          CALL section_vals_val_get(basis_section,"BASIS_SET",c_val=basis_name,error=error)
00518          gto_basis_section => section_vals_get_subs_vals(basis_section,"BASIS",error=error)
00519          CALL read_basis_set(ptable(zval)%symbol,basis,basis_name,basis_fn,&
00520                              gto_basis_section,error)
00521 
00522          ! initialize basis function on a radial grid
00523          nr = basis%grid%nr
00524          m  = MAXVAL(basis%nbas)
00525          ALLOCATE (basis%bf(nr,m,0:3),STAT=ierr)
00526          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00527          ALLOCATE (basis%dbf(nr,m,0:3),STAT=ierr)
00528          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00529          basis%bf =  0._dp
00530          basis%dbf = 0._dp
00531          DO l=0,3
00532            DO i=1,basis%nprim(l)
00533              al  = basis%am(i,l)
00534              DO k=1,nr
00535                rk  = basis%grid%rad(k)
00536                ear = EXP(-al*basis%grid%rad(k)**2)
00537                DO j=1,basis%nbas(l)
00538                  basis%bf(k,j,l) = basis%bf(k,j,l) + rk**l * ear*basis%cm(i,j,l)
00539                  basis%dbf(k,j,l) = basis%dbf(k,j,l) &
00540                    + ( REAL(l,dp)*rk**(l-1) - 2._dp*al*rk**(l+1) ) * ear*basis%cm(i,j,l)
00541                END DO
00542              END DO
00543            END DO
00544          END DO
00545       CASE (slater)
00546          basis%basis_type = STO_BASIS
00547          NULLIFY(num_slater)
00548          CALL section_vals_val_get(basis_section,"NUM_SLATER",i_vals=num_slater,error=error)
00549          IF ( num_slater(1) < 1 ) THEN
00550            CPAssert(.FALSE.,cp_failure_level,routineP,error,failure)
00551          ELSE
00552            basis%nbas = 0
00553            DO i=1,SIZE(num_slater)
00554              basis%nbas(i-1) = num_slater(i)
00555            END DO
00556            basis%nprim = basis%nbas
00557            m  = MAXVAL(basis%nbas)
00558            ALLOCATE (basis%as(m,0:3),basis%ns(m,0:3),STAT=ierr)
00559            CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00560            basis%as = 0._dp
00561            basis%ns = 0
00562            IF ( basis%nbas(0) > 0 ) THEN
00563              NULLIFY(expo)
00564              CALL section_vals_val_get(basis_section,"S_EXPONENTS",r_vals=expo,error=error)
00565              CPPostcondition(SIZE(expo)>=basis%nbas(0), cp_failure_level, routineP, error, failure)
00566              DO i=1,basis%nbas(0)
00567                basis%as(i,0) = expo(i)
00568              END DO
00569              NULLIFY(nqm)
00570              CALL section_vals_val_get(basis_section,"S_QUANTUM_NUMBERS",i_vals=nqm,error=error)
00571              CPPostcondition(SIZE(nqm)>=basis%nbas(0), cp_failure_level, routineP, error, failure)
00572              DO i=1,basis%nbas(0)
00573                basis%ns(i,0) = nqm(i)
00574              END DO
00575            END IF
00576            IF ( basis%nbas(1) > 0 ) THEN
00577              NULLIFY(expo)
00578              CALL section_vals_val_get(basis_section,"P_EXPONENTS",r_vals=expo,error=error)
00579              CPPostcondition(SIZE(expo)>=basis%nbas(1), cp_failure_level, routineP, error, failure)
00580              DO i=1,basis%nbas(1)
00581                basis%as(i,1) = expo(i)
00582              END DO
00583              NULLIFY(nqm)
00584              CALL section_vals_val_get(basis_section,"P_QUANTUM_NUMBERS",i_vals=nqm,error=error)
00585              CPPostcondition(SIZE(nqm)>=basis%nbas(1), cp_failure_level, routineP, error, failure)
00586              DO i=1,basis%nbas(1)
00587                basis%ns(i,1) = nqm(i)
00588              END DO
00589            END IF
00590            IF ( basis%nbas(2) > 0 ) THEN
00591              NULLIFY(expo)
00592              CALL section_vals_val_get(basis_section,"D_EXPONENTS",r_vals=expo,error=error)
00593              CPPostcondition(SIZE(expo)>=basis%nbas(2), cp_failure_level, routineP, error, failure)
00594              DO i=1,basis%nbas(2)
00595                basis%as(i,2) = expo(i)
00596              END DO
00597              NULLIFY(nqm)
00598              CALL section_vals_val_get(basis_section,"D_QUANTUM_NUMBERS",i_vals=nqm,error=error)
00599              CPPostcondition(SIZE(nqm)>=basis%nbas(2), cp_failure_level, routineP, error, failure)
00600              DO i=1,basis%nbas(2)
00601                basis%ns(i,2) = nqm(i)
00602              END DO
00603            END IF
00604            IF ( basis%nbas(3) > 0 ) THEN
00605              NULLIFY(expo)
00606              CALL section_vals_val_get(basis_section,"F_EXPONENTS",r_vals=expo,error=error)
00607              CPPostcondition(SIZE(expo)>=basis%nbas(3), cp_failure_level, routineP, error, failure)
00608              DO i=1,basis%nbas(3)
00609                basis%as(i,3) = expo(i)
00610              END DO
00611              NULLIFY(nqm)
00612              CALL section_vals_val_get(basis_section,"F_QUANTUM_NUMBERS",i_vals=nqm,error=error)
00613              CPPostcondition(SIZE(nqm)>=basis%nbas(3), cp_failure_level, routineP, error, failure)
00614              DO i=1,basis%nbas(3)
00615                basis%ns(i,3) = nqm(i)
00616              END DO
00617            END IF
00618          END IF
00619          ! initialize basis function on a radial grid
00620          nr = basis%grid%nr
00621          m  = MAXVAL(basis%nbas)
00622          ALLOCATE (basis%bf(nr,m,0:3),STAT=ierr)
00623          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00624          ALLOCATE (basis%dbf(nr,m,0:3),STAT=ierr)
00625          CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00626          basis%bf =  0._dp
00627          basis%dbf = 0._dp
00628          DO l=0,3
00629            DO i=1,basis%nbas(l)
00630              al  = basis%as(i,l)
00631              nl  = basis%ns(i,l)
00632              pf  = (2._dp*al)**nl * SQRT(2._dp*al/fac(2*nl))
00633              DO k=1,nr
00634                rk  = basis%grid%rad(k)
00635                ear = EXP(-al*rk)
00636                basis%bf(k,i,l) = pf * rk**(nl-1) * ear
00637                basis%dbf(k,i,l) = pf * ( REAL(nl-1,dp)/rk - al ) * rk**(nl-1) * ear
00638              END DO
00639            END DO
00640          END DO
00641       CASE (numerical)
00642          basis%basis_type = NUM_BASIS
00643          CPAssert(.FALSE.,cp_failure_level,routineP,error,failure)
00644     END SELECT
00645 
00646   END SUBROUTINE init_atom_basis
00647 
00648   SUBROUTINE release_atom_basis(basis,error)
00649     TYPE(atom_basis_type), INTENT(INOUT)     :: basis
00650     TYPE(cp_error_type), INTENT(inout)       :: error
00651 
00652     CHARACTER(len=*), PARAMETER :: routineN = 'release_atom_basis', 
00653       routineP = moduleN//':'//routineN
00654 
00655     INTEGER                                  :: stat
00656     LOGICAL                                  :: failure
00657 
00658     failure = .FALSE.
00659 
00660     IF(ASSOCIATED(basis%am)) THEN
00661        DEALLOCATE (basis%am,STAT=stat)
00662        CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00663     END IF
00664     IF(ASSOCIATED(basis%cm)) THEN
00665        DEALLOCATE (basis%cm,STAT=stat)
00666        CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00667     END IF
00668     IF(ASSOCIATED(basis%as)) THEN
00669        DEALLOCATE (basis%as,STAT=stat)
00670        CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00671     END IF
00672     IF(ASSOCIATED(basis%ns)) THEN
00673        DEALLOCATE (basis%ns,STAT=stat)
00674        CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00675     END IF
00676     IF(ASSOCIATED(basis%bf)) THEN
00677        DEALLOCATE (basis%bf,STAT=stat)
00678        CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00679     END IF
00680     IF(ASSOCIATED(basis%dbf)) THEN
00681        DEALLOCATE (basis%dbf,STAT=stat)
00682        CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00683     END IF
00684 
00685     CALL deallocate_grid_atom(basis%grid,error)
00686 
00687   END SUBROUTINE release_atom_basis
00688 ! *****************************************************************************
00689 
00690   SUBROUTINE create_atom_type(atom,error)
00691     TYPE(atom_type), POINTER                 :: atom
00692     TYPE(cp_error_type), INTENT(inout)       :: error
00693 
00694     CHARACTER(len=*), PARAMETER :: routineN = 'create_atom_type', 
00695       routineP = moduleN//':'//routineN
00696 
00697     INTEGER                                  :: ierr
00698     LOGICAL                                  :: failure
00699 
00700     failure = .FALSE.
00701 
00702     CPAssert(.NOT.ASSOCIATED(atom),cp_failure_level,routineP,error,failure)
00703 
00704     ALLOCATE(atom,STAT=ierr)
00705     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00706 
00707     NULLIFY(atom%xc_section)
00708 
00709   END SUBROUTINE create_atom_type
00710 
00711   SUBROUTINE release_atom_type(atom,error)
00712     TYPE(atom_type), POINTER                 :: atom
00713     TYPE(cp_error_type), INTENT(inout)       :: error
00714 
00715     CHARACTER(len=*), PARAMETER :: routineN = 'release_atom_type', 
00716       routineP = moduleN//':'//routineN
00717 
00718     INTEGER                                  :: ierr
00719     LOGICAL                                  :: failure
00720 
00721     failure = .FALSE.
00722 
00723     CPAssert(ASSOCIATED(atom),cp_failure_level,routineP,error,failure)
00724 
00725     NULLIFY(atom%basis)
00726     NULLIFY(atom%integrals)
00727     IF(ASSOCIATED(atom%state)) THEN
00728       DEALLOCATE(atom%state,STAT=ierr)
00729       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00730       NULLIFY(atom%state)
00731     END IF
00732     IF(ASSOCIATED(atom%orbitals)) THEN
00733       CALL release_atom_orbs(atom%orbitals,error)
00734     END IF
00735 
00736     DEALLOCATE(atom,STAT=ierr)
00737     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00738     NULLIFY(atom)
00739 
00740   END SUBROUTINE release_atom_type
00741 
00742   SUBROUTINE set_atom(atom,basis,state,integrals,orbitals,potential,zcore,pp_calc,&
00743     method_type,relativistic,coulomb_integral_type,exchange_integral_type,error)
00744     TYPE(atom_type), POINTER                 :: atom
00745     TYPE(atom_basis_type), OPTIONAL, POINTER :: basis
00746     TYPE(atom_state), OPTIONAL, POINTER      :: state
00747     TYPE(atom_integrals), OPTIONAL, POINTER  :: integrals
00748     TYPE(atom_orbitals), OPTIONAL, POINTER   :: orbitals
00749     TYPE(atom_potential_type), OPTIONAL, 
00750       POINTER                                :: potential
00751     INTEGER, INTENT(IN), OPTIONAL            :: zcore
00752     LOGICAL, INTENT(IN), OPTIONAL            :: pp_calc
00753     INTEGER, INTENT(IN), OPTIONAL            :: method_type, relativistic, 
00754                                                 coulomb_integral_type, 
00755                                                 exchange_integral_type
00756     TYPE(cp_error_type), INTENT(inout)       :: error
00757 
00758     CHARACTER(len=*), PARAMETER :: routineN = 'set_atom', 
00759       routineP = moduleN//':'//routineN
00760 
00761     LOGICAL                                  :: failure
00762 
00763     failure = .FALSE.
00764 
00765     CPAssert(ASSOCIATED(atom),cp_failure_level,routineP,error,failure)
00766 
00767     IF (.NOT.failure) THEN
00768       IF(PRESENT(basis)) atom%basis => basis
00769       IF(PRESENT(state)) atom%state => state
00770       IF(PRESENT(integrals)) atom%integrals => integrals
00771       IF(PRESENT(orbitals)) atom%orbitals => orbitals
00772       IF(PRESENT(potential)) atom%potential => potential
00773       IF(PRESENT(zcore)) atom%zcore = zcore
00774       IF(PRESENT(pp_calc)) atom%pp_calc = pp_calc
00775       IF(PRESENT(method_type)) atom%method_type = method_type
00776       IF(PRESENT(relativistic)) atom%relativistic = relativistic
00777       IF(PRESENT(coulomb_integral_type)) atom%coulomb_integral_type = coulomb_integral_type
00778       IF(PRESENT(exchange_integral_type)) atom%exchange_integral_type = exchange_integral_type
00779     END IF
00780 
00781   END SUBROUTINE set_atom
00782 
00783 ! *****************************************************************************
00784   SUBROUTINE create_atom_orbs(orbs,mbas,mo,error)
00785     TYPE(atom_orbitals), POINTER             :: orbs
00786     INTEGER, INTENT(IN)                      :: mbas, mo
00787     TYPE(cp_error_type), INTENT(inout)       :: error
00788 
00789     CHARACTER(len=*), PARAMETER :: routineN = 'create_atom_orbs', 
00790       routineP = moduleN//':'//routineN
00791 
00792     INTEGER                                  :: ierr
00793     LOGICAL                                  :: failure
00794 
00795     failure = .FALSE.
00796 
00797     CPAssert(.NOT.ASSOCIATED(orbs),cp_failure_level,routineP,error,failure)
00798 
00799     ALLOCATE(orbs,STAT=ierr)
00800     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00801 
00802     ALLOCATE(orbs%wfn(mbas,mo,0:3),orbs%wfna(mbas,mo,0:3),orbs%wfnb(mbas,mo,0:3),STAT=ierr)
00803     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00804     orbs%wfn = 0._dp
00805     orbs%wfna = 0._dp
00806     orbs%wfnb = 0._dp
00807 
00808     ALLOCATE(orbs%pmat(mbas,mbas,0:3),orbs%pmata(mbas,mbas,0:3),orbs%pmatb(mbas,mbas,0:3),STAT=ierr)
00809     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00810     orbs%pmat = 0._dp
00811     orbs%pmata = 0._dp
00812     orbs%pmatb = 0._dp
00813 
00814     ALLOCATE(orbs%ener(mo,0:3),orbs%enera(mo,0:3),orbs%enerb(mo,0:3),STAT=ierr)
00815     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00816     orbs%ener = 0._dp
00817     orbs%enera = 0._dp
00818     orbs%enerb = 0._dp
00819 
00820     ALLOCATE(orbs%refene(mo,0:3,2),orbs%refchg(mo,0:3,2),orbs%refnod(mo,0:3,2),STAT=ierr)
00821     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00822     orbs%refene = 0._dp
00823     orbs%refchg = 0._dp
00824     orbs%refnod = 0._dp
00825     ALLOCATE(orbs%wrefene(mo,0:3,2),orbs%wrefchg(mo,0:3,2),orbs%wrefnod(mo,0:3,2),STAT=ierr)
00826     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00827     orbs%wrefene = 0._dp
00828     orbs%wrefchg = 0._dp
00829     orbs%wrefnod = 0._dp
00830     ALLOCATE(orbs%crefene(mo,0:3,2),orbs%crefchg(mo,0:3,2),orbs%crefnod(mo,0:3,2),STAT=ierr)
00831     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00832     orbs%crefene = 0._dp
00833     orbs%crefchg = 0._dp
00834     orbs%crefnod = 0._dp
00835     ALLOCATE(orbs%rcmax(mo,0:3,2),STAT=ierr)
00836     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00837     orbs%rcmax = 0._dp
00838     ALLOCATE(orbs%wpsir0(mo,2),STAT=ierr)
00839     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00840     orbs%wpsir0 = 0._dp
00841 
00842   END SUBROUTINE create_atom_orbs
00843 
00844   SUBROUTINE release_atom_orbs(orbs,error)
00845     TYPE(atom_orbitals), POINTER             :: orbs
00846     TYPE(cp_error_type), INTENT(inout)       :: error
00847 
00848     CHARACTER(len=*), PARAMETER :: routineN = 'release_atom_orbs', 
00849       routineP = moduleN//':'//routineN
00850 
00851     INTEGER                                  :: ierr
00852     LOGICAL                                  :: failure
00853 
00854     failure = .FALSE.
00855 
00856     CPAssert(ASSOCIATED(orbs),cp_failure_level,routineP,error,failure)
00857 
00858     IF(ASSOCIATED(orbs%wfn)) THEN
00859       DEALLOCATE(orbs%wfn,orbs%wfna,orbs%wfnb,STAT=ierr)
00860       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00861       NULLIFY(orbs%wfn,orbs%wfna,orbs%wfnb)
00862     END IF
00863     IF(ASSOCIATED(orbs%pmat)) THEN
00864       DEALLOCATE(orbs%pmat,orbs%pmata,orbs%pmatb,STAT=ierr)
00865       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00866       NULLIFY(orbs%pmat,orbs%pmata,orbs%pmatb)
00867     END IF
00868     IF(ASSOCIATED(orbs%ener)) THEN
00869       DEALLOCATE(orbs%ener,orbs%enera,orbs%enerb,STAT=ierr)
00870       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00871       NULLIFY(orbs%ener,orbs%enera,orbs%enerb)
00872     END IF
00873     IF(ASSOCIATED(orbs%refene)) THEN
00874       DEALLOCATE(orbs%refene,STAT=ierr)
00875       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00876       NULLIFY(orbs%refene)
00877     END IF
00878     IF(ASSOCIATED(orbs%refchg)) THEN
00879       DEALLOCATE(orbs%refchg,STAT=ierr)
00880       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00881       NULLIFY(orbs%refchg)
00882     END IF
00883     IF(ASSOCIATED(orbs%refnod)) THEN
00884       DEALLOCATE(orbs%refnod,STAT=ierr)
00885       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00886       NULLIFY(orbs%refnod)
00887     END IF
00888     IF(ASSOCIATED(orbs%wrefene)) THEN
00889       DEALLOCATE(orbs%wrefene,STAT=ierr)
00890       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00891       NULLIFY(orbs%wrefene)
00892     END IF
00893     IF(ASSOCIATED(orbs%wrefchg)) THEN
00894       DEALLOCATE(orbs%wrefchg,STAT=ierr)
00895       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00896       NULLIFY(orbs%wrefchg)
00897     END IF
00898     IF(ASSOCIATED(orbs%wrefnod)) THEN
00899       DEALLOCATE(orbs%wrefnod,STAT=ierr)
00900       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00901       NULLIFY(orbs%wrefnod)
00902     END IF
00903     IF(ASSOCIATED(orbs%crefene)) THEN
00904       DEALLOCATE(orbs%crefene,STAT=ierr)
00905       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00906       NULLIFY(orbs%crefene)
00907     END IF
00908     IF(ASSOCIATED(orbs%crefchg)) THEN
00909       DEALLOCATE(orbs%crefchg,STAT=ierr)
00910       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00911       NULLIFY(orbs%crefchg)
00912     END IF
00913     IF(ASSOCIATED(orbs%crefnod)) THEN
00914       DEALLOCATE(orbs%crefnod,STAT=ierr)
00915       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00916       NULLIFY(orbs%crefnod)
00917     END IF
00918     IF(ASSOCIATED(orbs%rcmax)) THEN
00919       DEALLOCATE(orbs%rcmax,STAT=ierr)
00920       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00921       NULLIFY(orbs%rcmax)
00922     END IF
00923     IF(ASSOCIATED(orbs%wpsir0)) THEN
00924       DEALLOCATE(orbs%wpsir0,STAT=ierr)
00925       CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00926       NULLIFY(orbs%wpsir0)
00927     END IF
00928 
00929     DEALLOCATE(orbs,STAT=ierr)
00930     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00931     NULLIFY(orbs)
00932 
00933   END SUBROUTINE release_atom_orbs
00934 
00935 ! *****************************************************************************
00936   SUBROUTINE create_opmat(opmat,n,error)
00937     TYPE(opmat_type), POINTER                :: opmat
00938     INTEGER, DIMENSION(0:3), INTENT(IN)      :: n
00939     TYPE(cp_error_type), INTENT(inout)       :: error
00940 
00941     CHARACTER(len=*), PARAMETER :: routineN = 'create_opmat', 
00942       routineP = moduleN//':'//routineN
00943 
00944     INTEGER                                  :: ierr, m
00945     LOGICAL                                  :: failure = .FALSE.
00946 
00947     m=MAXVAL(n)
00948 
00949     CPPrecondition(.NOT.ASSOCIATED(opmat), cp_failure_level, routineP, error, failure)
00950 
00951     ALLOCATE(opmat,STAT=ierr)
00952     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00953 
00954     opmat%n = n
00955     ALLOCATE(opmat%op(m,m,0:3),STAT=ierr)
00956     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00957     opmat%op = 0._dp
00958 
00959   END SUBROUTINE create_opmat
00960 
00961   SUBROUTINE release_opmat(opmat,error)
00962     TYPE(opmat_type), POINTER                :: opmat
00963     TYPE(cp_error_type), INTENT(inout)       :: error
00964 
00965     CHARACTER(len=*), PARAMETER :: routineN = 'release_opmat', 
00966       routineP = moduleN//':'//routineN
00967 
00968     INTEGER                                  :: ierr
00969     LOGICAL                                  :: failure = .FALSE.
00970 
00971     CPPrecondition(ASSOCIATED(opmat), cp_failure_level, routineP, error, failure)
00972 
00973     opmat%n = 0
00974     DEALLOCATE(opmat%op,STAT=ierr)
00975     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00976 
00977     DEALLOCATE(opmat,STAT=ierr)
00978     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00979 
00980   END SUBROUTINE release_opmat
00981 
00982 ! *****************************************************************************
00983   SUBROUTINE create_opgrid(opgrid,grid,error)
00984     TYPE(opgrid_type), POINTER               :: opgrid
00985     TYPE(grid_atom_type), POINTER            :: grid
00986     TYPE(cp_error_type), INTENT(inout)       :: error
00987 
00988     CHARACTER(len=*), PARAMETER :: routineN = 'create_opgrid', 
00989       routineP = moduleN//':'//routineN
00990 
00991     INTEGER                                  :: ierr, nr
00992     LOGICAL                                  :: failure = .FALSE.
00993 
00994     CPPrecondition(.NOT.ASSOCIATED(opgrid), cp_failure_level, routineP, error, failure)
00995 
00996     ALLOCATE(opgrid,STAT=ierr)
00997     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
00998 
00999     opgrid%grid => grid
01000 
01001     nr = grid%nr
01002 
01003     ALLOCATE(opgrid%op(nr),STAT=ierr)
01004     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
01005     opgrid%op = 0._dp
01006 
01007   END SUBROUTINE create_opgrid
01008 
01009   SUBROUTINE release_opgrid(opgrid,error)
01010     TYPE(opgrid_type), POINTER               :: opgrid
01011     TYPE(cp_error_type), INTENT(inout)       :: error
01012 
01013     CHARACTER(len=*), PARAMETER :: routineN = 'release_opgrid', 
01014       routineP = moduleN//':'//routineN
01015 
01016     INTEGER                                  :: ierr
01017     LOGICAL                                  :: failure = .FALSE.
01018 
01019     CPPrecondition(ASSOCIATED(opgrid), cp_failure_level, routineP, error, failure)
01020 
01021     NULLIFY(opgrid%grid)
01022     DEALLOCATE(opgrid%op,STAT=ierr)
01023     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
01024 
01025     DEALLOCATE(opgrid,STAT=ierr)
01026     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
01027 
01028   END SUBROUTINE release_opgrid
01029 
01030 ! *****************************************************************************
01031   SUBROUTINE Clementi_geobas(zval,cval,aval,ngto,ival,error)
01032     INTEGER, INTENT(IN)                      :: zval
01033     REAL(dp), INTENT(OUT)                    :: cval, aval
01034     INTEGER, DIMENSION(0:3), INTENT(OUT)     :: ngto, ival
01035     TYPE(cp_error_type), INTENT(inout)       :: error
01036 
01037     CHARACTER(len=*), PARAMETER :: routineN = 'Clementi_geobas', 
01038       routineP = moduleN//':'//routineN
01039 
01040     LOGICAL                                  :: failure = .FALSE.
01041 
01042     ngto = 0
01043     ival = 0
01044     cval = 0._dp
01045     aval = 0._dp
01046 
01047     SELECT CASE (zval)
01048       CASE DEFAULT
01049         CPPrecondition(.FALSE., cp_failure_level, routineP, error, failure)
01050       CASE (1)   ! this is from the general geometrical basis and extended
01051         cval = 2.0_dp
01052         aval = 0.016_dp
01053         ngto(0) = 20
01054       CASE (2)
01055         cval = 2.14774520_dp
01056         aval = 0.04850670_dp
01057         ngto(0) = 20
01058       CASE (3)
01059         cval = 2.08932430_dp
01060         aval = 0.02031060_dp
01061         ngto(0) = 23
01062       CASE (4)
01063         cval = 2.09753060_dp
01064         aval = 0.03207070_dp
01065         ngto(0) = 23
01066       CASE (5)
01067         cval = 2.10343410_dp
01068         aval = 0.03591970_dp
01069         ngto(0) = 23
01070         ngto(1) = 16
01071       CASE (6)
01072         cval = 2.10662820_dp
01073         aval = 0.05292410_dp
01074         ngto(0) = 23
01075         ngto(1) = 16
01076       CASE (7)
01077         cval = 2.13743840_dp
01078         aval = 0.06291970_dp
01079         ngto(0) = 23
01080         ngto(1) = 16
01081       CASE (8)
01082         cval = 2.08687310_dp
01083         aval = 0.08350860_dp
01084         ngto(0) = 23
01085         ngto(1) = 16
01086       CASE (9)
01087         cval = 2.12318180_dp
01088         aval = 0.09899170_dp
01089         ngto(0) = 23
01090         ngto(1) = 16
01091       CASE (10)
01092         cval = 2.13164810_dp
01093         aval = 0.11485350_dp
01094         ngto(0) = 23
01095         ngto(1) = 16
01096       CASE (11)
01097         cval = 2.11413310_dp
01098         aval = 0.00922630_dp
01099         ngto(0) = 26
01100         ngto(1) = 16
01101         ival(1) = 4
01102       CASE (12)
01103         cval = 2.12183620_dp
01104         aval = 0.01215850_dp
01105         ngto(0) = 26
01106         ngto(1) = 16
01107         ival(1) = 4
01108       CASE (13)
01109         cval = 2.06073230_dp
01110         aval = 0.01449350_dp
01111         ngto(0) = 26
01112         ngto(1) = 20
01113         ival(0) = 1
01114       CASE (14)
01115         cval = 2.08563660_dp
01116         aval = 0.01861460_dp
01117         ngto(0) = 26
01118         ngto(1) = 20
01119         ival(0) = 1
01120       CASE (15)
01121         cval = 2.04879270_dp
01122         aval = 0.02147790_dp
01123         ngto(0) = 26
01124         ngto(1) = 20
01125         ival(0) = 1
01126       CASE (16)
01127         cval = 2.06216660_dp
01128         aval = 0.01978920_dp
01129         ngto(0) = 26
01130         ngto(1) = 20
01131         ival(0) = 1
01132       CASE (17)
01133         cval = 2.04628670_dp
01134         aval = 0.02451470_dp
01135         ngto(0) = 26
01136         ngto(1) = 20
01137         ival(0) = 1
01138       CASE (18)
01139         cval = 2.08675200_dp
01140         aval = 0.02635040_dp
01141         ngto(0) = 26
01142         ngto(1) = 20
01143         ival(0) = 1
01144       CASE (19)
01145         cval = 2.02715220_dp
01146         aval = 0.01822040_dp
01147         ngto(0) = 29
01148         ngto(1) = 20
01149         ival(1) = 2
01150       CASE (20)
01151         cval = 2.01465650_dp
01152         aval = 0.01646570_dp
01153         ngto(0) = 29
01154         ngto(1) = 20
01155         ival(1) = 2
01156       CASE (21)
01157         cval = 2.01605240_dp
01158         aval = 0.01254190_dp
01159         ngto(0) = 30
01160         ngto(1) = 20
01161         ngto(2) = 18
01162         ival(1) = 2
01163       CASE (22)
01164         cval = 2.01800000_dp
01165         aval = 0.01195490_dp
01166         ngto(0) = 30
01167         ngto(1) = 21
01168         ngto(2) = 17
01169         ival(1) = 2
01170         ival(2) = 1
01171       CASE (23)
01172         cval = 1.98803560_dp
01173         aval = 0.02492140_dp
01174         ngto(0) = 30
01175         ngto(1) = 21
01176         ngto(2) = 17
01177         ival(1) = 2
01178         ival(2) = 1
01179       CASE (24)
01180         cval = 1.98984000_dp
01181         aval = 0.02568400_dp
01182         ngto(0) = 30
01183         ngto(1) = 21
01184         ngto(2) = 17
01185         ival(1) = 2
01186         ival(2) = 1
01187       CASE (25)
01188         cval = 2.01694380_dp
01189         aval = 0.02664480_dp
01190         ngto(0) = 30
01191         ngto(1) = 21
01192         ngto(2) = 17
01193         ival(1) = 2
01194         ival(2) = 1
01195       CASE (26)
01196         cval = 2.01824090_dp
01197         aval = 0.01355000_dp
01198         ngto(0) = 30
01199         ngto(1) = 21
01200         ngto(2) = 17
01201         ival(1) = 2
01202         ival(2) = 1
01203       CASE (27)
01204         cval = 1.98359400_dp
01205         aval = 0.01702210_dp
01206         ngto(0) = 30
01207         ngto(1) = 21
01208         ngto(2) = 17
01209         ival(1) = 2
01210         ival(2) = 2
01211       CASE (28)
01212         cval = 1.96797340_dp
01213         aval = 0.02163180_dp
01214         ngto(0) = 30
01215         ngto(1) = 22
01216         ngto(2) = 17
01217         ival(1) = 3
01218         ival(2) = 2
01219       CASE (29)
01220         cval = 1.98955180_dp
01221         aval = 0.02304480_dp
01222         ngto(0) = 30
01223         ngto(1) = 20
01224         ngto(2) = 17
01225         ival(1) = 3
01226         ival(2) = 2
01227       CASE (30)
01228         cval = 1.98074320_dp
01229         aval = 0.02754320_dp
01230         ngto(0) = 30
01231         ngto(1) = 21
01232         ngto(2) = 17
01233         ival(1) = 3
01234         ival(2) = 2
01235       CASE (31)
01236         cval = 2.00551070_dp
01237         aval = 0.02005530_dp
01238         ngto(0) = 30
01239         ngto(1) = 23
01240         ngto(2) = 17
01241         ival(0) = 1
01242         ival(2) = 2
01243       CASE (32)
01244         cval = 2.00000030_dp
01245         aval = 0.02003000_dp
01246         ngto(0) = 30
01247         ngto(1) = 24
01248         ngto(2) = 17
01249         ival(0) = 1
01250         ival(2) = 2
01251       CASE (33)
01252         cval = 2.00609100_dp
01253         aval = 0.02055620_dp
01254         ngto(0) = 30
01255         ngto(1) = 23
01256         ngto(2) = 17
01257         ival(0) = 1
01258         ival(2) = 2
01259       CASE (34)
01260         cval = 2.00701000_dp
01261         aval = 0.02230400_dp
01262         ngto(0) = 30
01263         ngto(1) = 24
01264         ngto(2) = 17
01265         ival(0) = 1
01266         ival(2) = 2
01267       CASE (35)
01268         cval = 2.01508710_dp
01269         aval = 0.02685790_dp
01270         ngto(0) = 30
01271         ngto(1) = 24
01272         ngto(2) = 17
01273         ival(0) = 1
01274         ival(2) = 2
01275       CASE (36)
01276         cval = 2.01960430_dp
01277         aval = 0.02960430_dp
01278         ngto(0) = 30
01279         ngto(1) = 24
01280         ngto(2) = 17
01281         ival(0) = 1
01282         ival(2) = 2
01283       CASE (37)
01284         cval = 2.00031000_dp
01285         aval = 0.00768400_dp
01286         ngto(0) = 32
01287         ngto(1) = 25
01288         ngto(2) = 17
01289         ival(0) = 1
01290         ival(1) = 1
01291         ival(2) = 4
01292       CASE (38)
01293         cval = 1.99563960_dp
01294         aval = 0.01401940_dp
01295         ngto(0) = 33
01296         ngto(1) = 24
01297         ngto(2) = 17
01298         ival(1) = 1
01299         ival(2) = 4
01300       CASE (39)
01301         cval = 1.98971210_dp
01302         aval = 0.01558470_dp
01303         ngto(0) = 33
01304         ngto(1) = 24
01305         ngto(2) = 20
01306         ival(1) = 1
01307       CASE (40)
01308         cval = 1.97976190_dp
01309         aval = 0.01705520_dp
01310         ngto(0) = 33
01311         ngto(1) = 24
01312         ngto(2) = 20
01313         ival(1) = 1
01314       CASE (41)
01315         cval = 1.97989290_dp
01316         aval = 0.01527040_dp
01317         ngto(0) = 33
01318         ngto(1) = 24
01319         ngto(2) = 20
01320         ival(1) = 1
01321       CASE (42)
01322         cval = 1.97909240_dp
01323         aval = 0.01879720_dp
01324         ngto(0) = 32
01325         ngto(1) = 24
01326         ngto(2) = 20
01327         ival(1) = 1
01328       CASE (43)
01329         cval = 1.98508430_dp
01330         aval = 0.01497550_dp
01331         ngto(0) = 32
01332         ngto(1) = 24
01333         ngto(2) = 20
01334         ival(1) = 2
01335         ival(2) = 1
01336       CASE (44)
01337         cval = 1.98515010_dp
01338         aval = 0.01856670_dp
01339         ngto(0) = 32
01340         ngto(1) = 24
01341         ngto(2) = 20
01342         ival(1) = 2
01343         ival(2) = 1
01344       CASE (45)
01345         cval = 1.98502970_dp
01346         aval = 0.01487000_dp
01347         ngto(0) = 32
01348         ngto(1) = 24
01349         ngto(2) = 20
01350         ival(1) = 2
01351         ival(2) = 1
01352       CASE (46)
01353         cval = 1.97672850_dp
01354         aval = 0.01762500_dp
01355         ngto(0) = 30
01356         ngto(1) = 24
01357         ngto(2) = 20
01358         ival(0) = 2
01359         ival(1) = 2
01360         ival(2) = 1
01361       CASE (47)
01362         cval = 1.97862730_dp
01363         aval = 0.01863310_dp
01364         ngto(0) = 32
01365         ngto(1) = 24
01366         ngto(2) = 20
01367         ival(1) = 2
01368         ival(2) = 1
01369       CASE (48)
01370         cval = 1.97990020_dp
01371         aval = 0.01347150_dp
01372         ngto(0) = 33
01373         ngto(1) = 24
01374         ngto(2) = 20
01375         ival(1) = 2
01376         ival(2) = 2
01377       CASE (49)
01378         cval = 1.97979410_dp
01379         aval = 0.00890265_dp
01380         ngto(0) = 33
01381         ngto(1) = 27
01382         ngto(2) = 20
01383         ival(0) = 2
01384         ival(2) = 2
01385       CASE (50)
01386         cval = 1.98001000_dp
01387         aval = 0.00895215_dp
01388         ngto(0) = 33
01389         ngto(1) = 27
01390         ngto(2) = 20
01391         ival(0) = 2
01392         ival(2) = 2
01393       CASE (51)
01394         cval = 1.97979980_dp
01395         aval = 0.01490290_dp
01396         ngto(0) = 33
01397         ngto(1) = 26
01398         ngto(2) = 20
01399         ival(1) = 1
01400         ival(2) = 2
01401       CASE (52)
01402         cval = 1.98009310_dp
01403         aval = 0.01490390_dp
01404         ngto(0) = 33
01405         ngto(1) = 26
01406         ngto(2) = 20
01407         ival(1) = 1
01408         ival(2) = 2
01409       CASE (53)
01410         cval = 1.97794750_dp
01411         aval = 0.01425880_dp
01412         ngto(0) = 33
01413         ngto(1) = 26
01414         ngto(2) = 20
01415         ival(0) = 2
01416         ival(1) = 1
01417         ival(2) = 2
01418       CASE (54)
01419         cval = 1.97784450_dp
01420         aval = 0.01430130_dp
01421         ngto(0) = 33
01422         ngto(1) = 26
01423         ngto(2) = 20
01424         ival(0) = 2
01425         ival(1) = 1
01426         ival(2) = 2
01427       CASE (55)
01428         cval = 1.97784450_dp
01429         aval = 0.00499318_dp
01430         ngto(0) = 32
01431         ngto(1) = 25
01432         ngto(2) = 17
01433         ival(0) = 1
01434         ival(1) = 3
01435         ival(2) = 6
01436       CASE (56)
01437         cval = 1.97764820_dp
01438         aval = 0.00500392_dp
01439         ngto(0) = 32
01440         ngto(1) = 25
01441         ngto(2) = 17
01442         ival(0) = 1
01443         ival(1) = 3
01444         ival(2) = 6
01445       CASE (57)
01446         cval = 1.97765150_dp
01447         aval = 0.00557083_dp
01448         ngto(0) = 32
01449         ngto(1) = 25
01450         ngto(2) = 20
01451         ival(0) = 1
01452         ival(1) = 3
01453         ival(2) = 3
01454       CASE (58)
01455         cval = 1.97768750_dp
01456         aval = 0.00547531_dp
01457         ngto(0) = 32
01458         ngto(1) = 25
01459         ngto(2) = 20
01460         ngto(3) = 16
01461         ival(0) = 1
01462         ival(1) = 3
01463         ival(2) = 3
01464         ival(3) = 3
01465       CASE (59)
01466         cval = 1.96986600_dp
01467         aval = 0.00813143_dp
01468         ngto(0) = 32
01469         ngto(1) = 25
01470         ngto(2) = 17
01471         ngto(3) = 16
01472         ival(0) = 1
01473         ival(1) = 3
01474         ival(2) = 6
01475         ival(3) = 4
01476       CASE (60)
01477         cval = 1.97765720_dp
01478         aval = 0.00489201_dp
01479         ngto(0) = 32
01480         ngto(1) = 25
01481         ngto(2) = 17
01482         ngto(3) = 16
01483         ival(0) = 1
01484         ival(1) = 3
01485         ival(2) = 6
01486         ival(3) = 4
01487       CASE (61)
01488         cval = 1.97768120_dp
01489         aval = 0.00499000_dp
01490         ngto(0) = 32
01491         ngto(1) = 25
01492         ngto(2) = 17
01493         ngto(3) = 16
01494         ival(0) = 1
01495         ival(1) = 3
01496         ival(2) = 6
01497         ival(3) = 4
01498       CASE (62)
01499         cval = 1.97745700_dp
01500         aval = 0.00615587_dp
01501         ngto(0) = 32
01502         ngto(1) = 25
01503         ngto(2) = 17
01504         ngto(3) = 16
01505         ival(0) = 1
01506         ival(1) = 3
01507         ival(2) = 6
01508         ival(3) = 4
01509       CASE (63)
01510         cval = 1.97570240_dp
01511         aval = 0.00769959_dp
01512         ngto(0) = 32
01513         ngto(1) = 25
01514         ngto(2) = 17
01515         ngto(3) = 16
01516         ival(0) = 1
01517         ival(1) = 3
01518         ival(2) = 6
01519         ival(3) = 4
01520       CASE (64)
01521         cval = 1.97629350_dp
01522         aval = 0.00706610_dp
01523         ngto(0) = 32
01524         ngto(1) = 25
01525         ngto(2) = 20
01526         ngto(3) = 16
01527         ival(0) = 1
01528         ival(1) = 3
01529         ival(2) = 3
01530         ival(3) = 4
01531       CASE (65)
01532         cval = 1.96900000_dp
01533         aval = 0.01019150_dp
01534         ngto(0) = 32
01535         ngto(1) = 26
01536         ngto(2) = 18
01537         ngto(3) = 16
01538         ival(0) = 1
01539         ival(1) = 3
01540         ival(2) = 6
01541         ival(3) = 4
01542       CASE (66)
01543         cval = 1.97350000_dp
01544         aval = 0.01334320_dp
01545         ngto(0) = 33
01546         ngto(1) = 26
01547         ngto(2) = 18
01548         ngto(3) = 16
01549         ival(0) = 1
01550         ival(1) = 3
01551         ival(2) = 6
01552         ival(3) = 4
01553       CASE (67)
01554         cval = 1.97493000_dp
01555         aval = 0.01331360_dp
01556         ngto(0) = 32
01557         ngto(1) = 24
01558         ngto(2) = 17
01559         ngto(3) = 14
01560         ival(1) = 2
01561         ival(2) = 5
01562         ival(3) = 4
01563       CASE (68)
01564         cval = 1.97597670_dp
01565         aval = 0.01434040_dp
01566         ngto(0) = 32
01567         ngto(1) = 24
01568         ngto(2) = 17
01569         ngto(3) = 14
01570         ival(0) = 0
01571         ival(1) = 2
01572         ival(2) = 5
01573         ival(3) = 4
01574       CASE (69)
01575         cval = 1.97809240_dp
01576         aval = 0.01529430_dp
01577         ngto(0) = 32
01578         ngto(1) = 24
01579         ngto(2) = 17
01580         ngto(3) = 14
01581         ival(0) = 0
01582         ival(1) = 2
01583         ival(2) = 5
01584         ival(3) = 4
01585       CASE (70)
01586         cval = 1.97644360_dp
01587         aval = 0.01312770_dp
01588         ngto(0) = 32
01589         ngto(1) = 24
01590         ngto(2) = 17
01591         ngto(3) = 14
01592         ival(0) = 0
01593         ival(1) = 2
01594         ival(2) = 5
01595         ival(3) = 4
01596       CASE (71)
01597         cval = 1.96998000_dp
01598         aval = 0.01745150_dp
01599         ngto(0) = 31
01600         ngto(1) = 24
01601         ngto(2) = 20
01602         ngto(3) = 14
01603         ival(0) = 1
01604         ival(1) = 2
01605         ival(2) = 2
01606         ival(3) = 4
01607       CASE (72)
01608         cval = 1.97223830_dp
01609         aval = 0.01639750_dp
01610         ngto(0) = 31
01611         ngto(1) = 24
01612         ngto(2) = 20
01613         ngto(3) = 14
01614         ival(0) = 1
01615         ival(1) = 2
01616         ival(2) = 2
01617         ival(3) = 4
01618       CASE (73)
01619         cval = 1.97462110_dp
01620         aval = 0.01603680_dp
01621         ngto(0) = 31
01622         ngto(1) = 24
01623         ngto(2) = 20
01624         ngto(3) = 14
01625         ival(0) = 1
01626         ival(1) = 2
01627         ival(2) = 2
01628         ival(3) = 4
01629       CASE (74)
01630         cval = 1.97756000_dp
01631         aval = 0.02030570_dp
01632         ngto(0) = 31
01633         ngto(1) = 24
01634         ngto(2) = 20
01635         ngto(3) = 14
01636         ival(0) = 1
01637         ival(1) = 2
01638         ival(2) = 2
01639         ival(3) = 4
01640       CASE (75)
01641         cval = 1.97645760_dp
01642         aval = 0.02057180_dp
01643         ngto(0) = 31
01644         ngto(1) = 24
01645         ngto(2) = 20
01646         ngto(3) = 14
01647         ival(0) = 1
01648         ival(1) = 2
01649         ival(2) = 2
01650         ival(3) = 4
01651       CASE (76)
01652         cval = 1.97725820_dp
01653         aval = 0.02058210_dp
01654         ngto(0) = 32
01655         ngto(1) = 24
01656         ngto(2) = 20
01657         ngto(3) = 15
01658         ival(0) = 0
01659         ival(1) = 2
01660         ival(2) = 2
01661         ival(3) = 4
01662       CASE (77)
01663         cval = 1.97749380_dp
01664         aval = 0.02219380_dp
01665         ngto(0) = 32
01666         ngto(1) = 24
01667         ngto(2) = 20
01668         ngto(3) = 15
01669         ival(0) = 0
01670         ival(1) = 2
01671         ival(2) = 2
01672         ival(3) = 4
01673       CASE (78)
01674         cval = 1.97946280_dp
01675         aval = 0.02216280_dp
01676         ngto(0) = 32
01677         ngto(1) = 24
01678         ngto(2) = 20
01679         ngto(3) = 15
01680         ival(0) = 0
01681         ival(1) = 2
01682         ival(2) = 2
01683         ival(3) = 4
01684       CASE (79)
01685         cval = 1.97852130_dp
01686         aval = 0.02168500_dp
01687         ngto(0) = 32
01688         ngto(1) = 24
01689         ngto(2) = 20
01690         ngto(3) = 15
01691         ival(0) = 0
01692         ival(1) = 2
01693         ival(2) = 2
01694         ival(3) = 4
01695       CASE (80)
01696         cval = 1.98045190_dp
01697         aval = 0.02177860_dp
01698         ngto(0) = 32
01699         ngto(1) = 24
01700         ngto(2) = 20
01701         ngto(3) = 15
01702         ival(0) = 0
01703         ival(1) = 2
01704         ival(2) = 2
01705         ival(3) = 4
01706       CASE (81)
01707         cval = 1.97000000_dp
01708         aval = 0.02275000_dp
01709         ngto(0) = 31
01710         ngto(1) = 25
01711         ngto(2) = 18
01712         ngto(3) = 13
01713         ival(0) = 1
01714         ival(1) = 0
01715         ival(2) = 3
01716         ival(3) = 6
01717       CASE (82)
01718         cval = 1.97713580_dp
01719         aval = 0.02317030_dp
01720         ngto(0) = 31
01721         ngto(1) = 27
01722         ngto(2) = 18
01723         ngto(3) = 13
01724         ival(0) = 1
01725         ival(1) = 0
01726         ival(2) = 3
01727         ival(3) = 6
01728       CASE (83)
01729         cval = 1.97537880_dp
01730         aval = 0.02672860_dp
01731         ngto(0) = 32
01732         ngto(1) = 27
01733         ngto(2) = 17
01734         ngto(3) = 13
01735         ival(0) = 1
01736         ival(1) = 0
01737         ival(2) = 3
01738         ival(3) = 6
01739       CASE (84)
01740         cval = 1.97545360_dp
01741         aval = 0.02745360_dp
01742         ngto(0) = 31
01743         ngto(1) = 27
01744         ngto(2) = 17
01745         ngto(3) = 13
01746         ival(0) = 1
01747         ival(1) = 0
01748         ival(2) = 3
01749         ival(3) = 6
01750       CASE (85)
01751         cval = 1.97338370_dp
01752         aval = 0.02616310_dp
01753         ngto(0) = 32
01754         ngto(1) = 27
01755         ngto(2) = 19
01756         ngto(3) = 13
01757         ival(0) = 1
01758         ival(1) = 0
01759         ival(2) = 3
01760         ival(3) = 6
01761       CASE (86)
01762         cval = 1.97294240_dp
01763         aval = 0.02429220_dp
01764         ngto(0) = 32
01765         ngto(1) = 27
01766         ngto(2) = 19
01767         ngto(3) = 13
01768         ival(0) = 1
01769         ival(1) = 0
01770         ival(2) = 3
01771         ival(3) = 6
01772       CASE (87:106) ! these numbers are a educated guess
01773         cval = 1.98000000_dp
01774         aval = 0.01400000_dp
01775         ngto(0) = 34
01776         ngto(1) = 28
01777         ngto(2) = 20
01778         ngto(3) = 15
01779         ival(0) = 0
01780         ival(1) = 0
01781         ival(2) = 3
01782         ival(3) = 6
01783     END SELECT
01784 
01785   END SUBROUTINE Clementi_geobas
01786 ! *****************************************************************************
01787   SUBROUTINE read_basis_set(element_symbol,basis,basis_set_name,basis_set_file,&
01788                             basis_section,error)
01789 
01790     CHARACTER(LEN=*), INTENT(IN)             :: element_symbol
01791     TYPE(atom_basis_type), INTENT(INOUT)     :: basis
01792     CHARACTER(LEN=*), INTENT(IN)             :: basis_set_name, basis_set_file
01793     TYPE(section_vals_type), POINTER         :: basis_section
01794     TYPE(cp_error_type), INTENT(inout)       :: error
01795 
01796     CHARACTER(len=*), PARAMETER :: routineN = 'read_basis_set', 
01797       routineP = moduleN//':'//routineN
01798     INTEGER, PARAMETER                       :: maxpri = 40, maxset = 20
01799 
01800     CHARACTER(len=20*default_string_length)  :: line_att
01801     CHARACTER(LEN=240)                       :: line
01802     CHARACTER(LEN=242)                       :: line2
01803     CHARACTER(LEN=LEN(basis_set_name))       :: bsname
01804     CHARACTER(LEN=LEN(basis_set_name)+2)     :: bsname2
01805     CHARACTER(LEN=LEN(element_symbol))       :: symbol
01806     CHARACTER(LEN=LEN(element_symbol)+2)     :: symbol2
01807     INTEGER                                  :: i, ierr, ii, ipgf, irep, 
01808                                                 iset, ishell, j, k, lshell, 
01809                                                 nj, nmin, ns, nset, strlen1, 
01810                                                 strlen2
01811     INTEGER, DIMENSION(maxpri, maxset)       :: l
01812     INTEGER, DIMENSION(maxset)               :: lmax, lmin, n, npgf, nshell
01813     LOGICAL                                  :: failure, found, is_ok, match, 
01814                                                 read_from_input
01815     REAL(dp)                                 :: expzet, gcca, prefac, zeta
01816     REAL(dp), 
01817       DIMENSION(maxpri, maxpri, maxset)      :: gcc
01818     REAL(dp), DIMENSION(maxpri, maxset)      :: zet
01819     TYPE(cp_parser_type), POINTER            :: parser
01820     TYPE(cp_sll_val_type), POINTER           :: list
01821     TYPE(val_type), POINTER                  :: val
01822 
01823     failure = .FALSE.
01824 
01825     bsname = basis_set_name
01826     symbol = element_symbol
01827     irep   = 0
01828 
01829     nset = 0
01830     lmin = 0
01831     lmax = 0
01832     npgf = 0
01833     n = 0
01834     l = 0
01835     zet = 0._dp
01836     gcc = 0._dp
01837 
01838     read_from_input = .FALSE.
01839     CALL section_vals_get(basis_section,explicit=read_from_input, error=error)
01840     IF (read_from_input) THEN
01841       NULLIFY(list,val)
01842       CALL section_vals_list_get(basis_section,"_DEFAULT_KEYWORD_",list=list,error=error)
01843       CALL uppercase(symbol)
01844       CALL uppercase(bsname)
01845       is_ok=cp_sll_val_next(list,val,error=error)
01846       CPPrecondition(is_ok, cp_failure_level, routineP, error, failure)
01847       CALL val_get(val,c_val=line_att,error=error)
01848       READ(line_att,*) nset
01849       CPPrecondition(nset <= maxset, cp_failure_level, routineP, error, failure)
01850       DO iset=1,nset
01851         is_ok=cp_sll_val_next(list,val,error=error)
01852         CPPrecondition(is_ok, cp_failure_level, routineP, error, failure)
01853         CALL val_get(val,c_val=line_att,error=error)
01854         READ(line_att,*) n(iset)
01855         CALL remove_word(line_att)
01856         READ(line_att,*) lmin(iset)
01857         CALL remove_word(line_att)
01858         READ(line_att,*) lmax(iset)
01859         CALL remove_word(line_att)
01860         READ(line_att,*) npgf(iset)
01861         CALL remove_word(line_att)
01862         CPPrecondition(npgf(iset) <= maxpri, cp_failure_level, routineP, error, failure)
01863         nshell(iset) = 0
01864         DO lshell=lmin(iset),lmax(iset)
01865           nmin = n(iset) + lshell - lmin(iset)
01866           READ(line_att,*) ishell
01867           CALL remove_word(line_att)
01868           nshell(iset) = nshell(iset) + ishell
01869           DO i=1,ishell
01870             l(nshell(iset)-ishell+i,iset) = lshell
01871           END DO
01872         END DO
01873         CPPrecondition(LEN_TRIM(line_att)==0, cp_failure_level, routineP, error, failure)
01874         DO ipgf=1,npgf(iset)
01875           is_ok=cp_sll_val_next(list,val,error=error)
01876           CPPrecondition(is_ok, cp_failure_level, routineP, error, failure)
01877           CALL val_get(val,c_val=line_att,error=error)
01878           READ(line_att,*) zet(ipgf,iset),(gcc(ipgf,ishell,iset),ishell=1,nshell(iset))
01879         END DO
01880       END DO
01881     ELSE
01882       NULLIFY(parser)
01883       CALL parser_create(parser,basis_set_file,error=error)
01884       ! Search for the requested basis set in the basis set file
01885       ! until the basis set is found or the end of file is reached
01886       search_loop: DO
01887         CALL parser_search_string(parser,TRIM(bsname),.TRUE.,found,line,error=error)
01888         IF (found) THEN
01889           CALL uppercase(symbol)
01890           CALL uppercase(bsname)
01891           match = .FALSE.
01892           CALL uppercase(line)
01893           ! Check both the element symbol and the basis set name
01894           line2 = " "//line//" "
01895           symbol2 = " "//TRIM(symbol)//" "
01896           bsname2 = " "//TRIM(bsname)//" "
01897           strlen1 = LEN_TRIM(symbol2) + 1
01898           strlen2 = LEN_TRIM(bsname2) + 1
01899 
01900           IF ( (INDEX(line2,symbol2(:strlen1)) > 0).AND.&
01901                (INDEX(line2,bsname2(:strlen2)) > 0) ) match = .TRUE.
01902 
01903           IF (match) THEN
01904             ! Read the basis set information
01905             CALL parser_get_object(parser,nset,newline=.TRUE.,error=error)
01906             CPPrecondition(nset <= maxset, cp_failure_level, routineP, error, failure)
01907             DO iset=1,nset
01908               CALL parser_get_object(parser,n(iset),newline=.TRUE.,error=error)
01909               CALL parser_get_object(parser,lmin(iset),error=error)
01910               CALL parser_get_object(parser,lmax(iset),error=error)
01911               CALL parser_get_object(parser,npgf(iset),error=error)
01912               CPPrecondition(npgf(iset) <= maxpri, cp_failure_level, routineP, error, failure)
01913               nshell(iset) = 0
01914               DO lshell=lmin(iset),lmax(iset)
01915                 nmin = n(iset) + lshell - lmin(iset)
01916                 CALL parser_get_object(parser,ishell,error=error)
01917                 nshell(iset) = nshell(iset) + ishell
01918                 DO i=1,ishell
01919                   l(nshell(iset)-ishell+i,iset) = lshell
01920                 END DO
01921               END DO
01922               DO ipgf=1,npgf(iset)
01923                 CALL parser_get_object(parser,zet(ipgf,iset),newline=.TRUE.,error=error)
01924                 DO ishell=1,nshell(iset)
01925                   CALL parser_get_object(parser,gcc(ipgf,ishell,iset),error=error)
01926                 END DO
01927               END DO
01928             END DO
01929 
01930             EXIT search_loop
01931 
01932           END IF
01933         ELSE
01934           ! Stop program, if the end of file is reached
01935           CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure)
01936         END IF
01937 
01938       END DO search_loop
01939 
01940       CALL parser_release(parser,error=error)
01941     END IF
01942 
01943     ! fill in the basis data structures
01944     basis%nprim = 0
01945     basis%nbas = 0
01946     DO i=1,nset
01947       DO j=lmin(i),lmax(i)
01948         basis%nprim(j) = basis%nprim(j) + npgf(i)
01949       END DO
01950       DO j=1,nshell(i)
01951         k = l(j,i)
01952         basis%nbas(k) = basis%nbas(k)+1
01953       END DO
01954     END DO
01955 
01956     nj = MAXVAL(basis%nprim)
01957     ns = MAXVAL(basis%nbas)
01958     ALLOCATE (basis%am(nj,0:3),STAT=ierr)
01959     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
01960     basis%am = 0._dp
01961     ALLOCATE (basis%cm(nj,ns,0:3),STAT=ierr)
01962     CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure)
01963     basis%cm = 0._dp
01964 
01965     DO j=0,3
01966       nj = 0
01967       ns = 0
01968       DO i=1,nset
01969         IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
01970           DO ipgf = 1,npgf(i)
01971             basis%am(nj+ipgf,j) = zet(ipgf,i)
01972           END DO
01973           DO ii=1,nshell(i)
01974             IF ( l(ii,i) == j ) THEN
01975               ns = ns + 1
01976               DO ipgf=1,npgf(i)
01977                 basis%cm(nj+ipgf,ns,j) = gcc(ipgf,ii,i)
01978               END DO
01979             END IF
01980           END DO
01981           nj = nj + npgf(i)
01982         END IF
01983       END DO
01984     END DO
01985 
01986     ! Normalization
01987     DO j=0,3
01988       expzet = 0.25_dp*REAL(2*j + 3,dp)
01989       prefac = SQRT( SQRT(pi)/2._dp**(j+2)*dfac(2*j+1) )
01990       DO ipgf=1,basis%nprim(j)
01991         DO ii=1,basis%nbas(j)
01992           gcca = basis%cm(ipgf,ii,j)
01993           zeta = 2._dp*basis%am(ipgf,j)
01994           basis%cm(ipgf,ii,j) = zeta**expzet*gcca/prefac
01995         END DO
01996       END DO
01997     END DO
01998 
01999   END SUBROUTINE read_basis_set
02000 
02001 ! *****************************************************************************
02002   SUBROUTINE read_atom_opt_section(optimization,opt_section,error)
02003     TYPE(atom_optimization_type), 
02004       INTENT(INOUT)                          :: optimization
02005     TYPE(section_vals_type), POINTER         :: opt_section
02006     TYPE(cp_error_type), INTENT(inout)       :: error
02007 
02008     CHARACTER(len=*), PARAMETER :: routineN = 'read_atom_opt_section', 
02009       routineP = moduleN//':'//routineN
02010 
02011     INTEGER                                  :: miter, ndiis
02012     REAL(KIND=dp)                            :: damp, eps_diis, eps_scf
02013 
02014     CALL section_vals_val_get(opt_section,"MAX_ITER",i_val=miter,error=error)
02015     CALL section_vals_val_get(opt_section,"EPS_SCF",r_val=eps_scf,error=error)
02016     CALL section_vals_val_get(opt_section,"N_DIIS",i_val=ndiis,error=error)
02017     CALL section_vals_val_get(opt_section,"EPS_DIIS",r_val=eps_diis,error=error)
02018     CALL section_vals_val_get(opt_section,"DAMPING",r_val=damp,error=error)
02019 
02020     optimization%max_iter = miter
02021     optimization%eps_scf  = eps_scf
02022     optimization%n_diis   = ndiis
02023     optimization%eps_diis = eps_diis
02024     optimization%damping  = damp
02025 
02026   END SUBROUTINE read_atom_opt_section
02027 ! *****************************************************************************
02028   SUBROUTINE init_atom_potential(potential,potential_section,zval,error)
02029     TYPE(atom_potential_type), INTENT(INOUT) :: potential
02030     TYPE(section_vals_type), POINTER         :: potential_section
02031     INTEGER, INTENT(IN)                      :: zval
02032     TYPE(cp_error_type), INTENT(inout)       :: error
02033 
02034     CHARACTER(len=*), PARAMETER :: routineN = 'init_atom_potential', 
02035       routineP = moduleN//':'//routineN
02036 
02037     CHARACTER(LEN=default_string_length)     :: pseudo_fn, pseudo_name
02038     LOGICAL                                  :: failure
02039     REAL(dp), DIMENSION(:), POINTER          :: convals
02040     TYPE(section_vals_type), POINTER         :: gth_potential_section
02041 
02042     failure = .FALSE.
02043 
02044     IF ( zval > 0 ) THEN
02045       CALL section_vals_val_get(potential_section,"PSEUDO_TYPE",i_val=potential%ppot_type,error=error)
02046 
02047       SELECT CASE (potential%ppot_type)
02048         CASE (gth_pseudo)
02049            CALL section_vals_val_get(potential_section,"POTENTIAL_FILE_NAME",c_val=pseudo_fn,error=error)
02050            CALL section_vals_val_get(potential_section,"POTENTIAL_NAME",c_val=pseudo_name,error=error)
02051            gth_potential_section => section_vals_get_subs_vals(potential_section,"GTH_POTENTIAL",error=error)
02052            CALL read_gth_potential(ptable(zval)%symbol,potential%gth_pot,&
02053                  pseudo_name,pseudo_fn,gth_potential_section,error)
02054         CASE (no_pseudo)
02055           ! do nothing
02056         CASE DEFAULT
02057           CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure)
02058       END SELECT
02059     ELSE
02060       potential%ppot_type = no_pseudo
02061     END IF
02062 
02063     ! confinement
02064     NULLIFY(convals)
02065     CALL section_vals_val_get(potential_section,"CONFINEMENT",r_vals=convals,error=error)
02066     IF ( SIZE (convals) >= 1 ) THEN
02067       IF ( convals(1) > 0._dp ) THEN
02068         potential%confinement = .TRUE.
02069         potential%acon = convals(1)
02070         IF ( SIZE (convals) >= 2 ) THEN
02071           potential%rcon = convals(2)
02072         ELSE
02073           potential%rcon = 4._dp
02074         END IF
02075         IF ( SIZE (convals) >= 3 ) THEN
02076           potential%ncon = NINT(convals(3))
02077         ELSE
02078           potential%ncon = 2
02079         END IF
02080       ELSE
02081         potential%confinement = .FALSE.
02082       END IF
02083     ELSE
02084       potential%confinement = .FALSE.
02085     END IF
02086 
02087   END SUBROUTINE init_atom_potential
02088 ! *****************************************************************************
02089   SUBROUTINE release_atom_potential(potential,error)
02090     TYPE(atom_potential_type), INTENT(INOUT) :: potential
02091     TYPE(cp_error_type), INTENT(inout)       :: error
02092 
02093     CHARACTER(len=*), PARAMETER :: routineN = 'release_atom_potential', 
02094       routineP = moduleN//':'//routineN
02095 
02096     potential%confinement = .FALSE.
02097 
02098   END SUBROUTINE release_atom_potential
02099 ! *****************************************************************************
02100   SUBROUTINE read_gth_potential(element_symbol,potential,pseudo_name,pseudo_file,&
02101                             potential_section,error)
02102 
02103     CHARACTER(LEN=*), INTENT(IN)             :: element_symbol
02104     TYPE(atom_gthpot_type), INTENT(INOUT)    :: potential
02105     CHARACTER(LEN=*), INTENT(IN)             :: pseudo_name, pseudo_file
02106     TYPE(section_vals_type), POINTER         :: potential_section
02107     TYPE(cp_error_type), INTENT(inout)       :: error
02108 
02109     CHARACTER(len=*), PARAMETER :: routineN = 'read_gth_potential', 
02110       routineP = moduleN//':'//routineN
02111 
02112     CHARACTER(LEN=240)                       :: line
02113     CHARACTER(LEN=242)                       :: line2
02114     CHARACTER(len=5*default_string_length)   :: line_att
02115     CHARACTER(LEN=LEN(element_symbol))       :: symbol
02116     CHARACTER(LEN=LEN(element_symbol)+2)     :: symbol2
02117     CHARACTER(LEN=LEN(pseudo_name))          :: apname
02118     CHARACTER(LEN=LEN(pseudo_name)+2)        :: apname2
02119     INTEGER                                  :: i, ic, ipot, j, l, nlmax, 
02120                                                 strlen1, strlen2
02121     INTEGER, DIMENSION(0:4)                  :: elec_conf
02122     LOGICAL                                  :: failure, found, is_ok, match, 
02123                                                 read_from_input
02124     TYPE(cp_parser_type), POINTER            :: parser
02125     TYPE(cp_sll_val_type), POINTER           :: list
02126     TYPE(val_type), POINTER                  :: val
02127 
02128     failure = .FALSE.
02129 
02130     elec_conf = 0
02131 
02132     apname = pseudo_name
02133     symbol = element_symbol
02134 
02135     potential%symbol = symbol
02136     potential%pname = apname
02137     potential%econf = 0
02138     potential%rc = 0._dp
02139     potential%ncl = 0
02140     potential%cl = 0._dp
02141     potential%nl = 0
02142     potential%rcnl = 0._dp
02143     potential%hnl = 0._dp
02144 
02145     potential%lpotextended = .FALSE.
02146     potential%lsdpot = .FALSE.
02147     potential%nlcc = .FALSE.
02148     potential%nexp_lpot = 0
02149     potential%nexp_lsd = 0
02150     potential%nexp_nlcc = 0
02151 
02152     read_from_input = .FALSE.
02153     CALL section_vals_get(potential_section,explicit=read_from_input, error=error)
02154     IF (read_from_input) THEN
02155       CALL section_vals_list_get(potential_section,"_DEFAULT_KEYWORD_",list=list,error=error)
02156       CALL uppercase(symbol)
02157       CALL uppercase(apname)
02158       ! Read the electronic configuration, not used here
02159       l = 0
02160       is_ok=cp_sll_val_next(list,val,error=error)
02161       CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02162       CALL val_get(val,c_val=line_att,error=error)
02163       READ(line_att,*) elec_conf(l)
02164       CALL remove_word(line_att)
02165       DO WHILE (LEN_TRIM(line_att) /= 0)
02166         l = l + 1
02167         READ(line_att,*)elec_conf(l)
02168         CALL remove_word(line_att)
02169       END DO
02170       potential%econf(0:3) = elec_conf(0:3)
02171       potential%zion = REAL ( SUM(elec_conf), dp )
02172       ! Read r(loc) to define the exponent of the core charge
02173       is_ok=cp_sll_val_next(list,val,error=error)
02174       CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02175       CALL val_get(val,c_val=line_att,error=error)
02176       READ(line_att,*) potential%rc
02177       CALL remove_word(line_att)
02178       ! Read the parameters for the local part of the GTH pseudopotential (ppl)
02179       READ(line_att,*) potential%ncl
02180       CALL remove_word(line_att)
02181       DO i=1,potential%ncl
02182         READ(line_att,*) potential%cl(i)
02183         CALL remove_word(line_att)
02184       END DO
02185       ! Check for the next entry: LPOT, NLCC, LSD, or ppnl
02186       DO
02187         is_ok=cp_sll_val_next(list,val,error=error)
02188         CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02189         CALL val_get(val,c_val=line_att,error=error)
02190         IF(INDEX(line_att,"LPOT") /= 0) THEN
02191           potential%lpotextended = .TRUE.
02192           CALL remove_word(line_att)
02193           READ(line_att,*) potential%nexp_lpot
02194           DO ipot=1,potential%nexp_lpot
02195              is_ok=cp_sll_val_next(list,val,error=error)
02196              CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02197              CALL val_get(val,c_val=line_att,error=error)
02198              READ(line_att,*) potential%alpha_lpot(ipot)
02199              CALL remove_word(line_att)
02200              READ(line_att,*) potential%nct_lpot(ipot)
02201              CALL remove_word(line_att)
02202              DO ic=1,potential%nct_lpot(ipot)
02203                READ(line_att,*) potential%cval_lpot(ic,ipot)
02204                CALL remove_word(line_att)
02205              END DO
02206           END DO
02207         ELSEIF(INDEX(line_att,"NLCC") /= 0) THEN
02208           potential%nlcc = .TRUE.
02209           CALL remove_word(line_att)
02210           READ(line_att,*) potential%nexp_nlcc
02211           DO ipot=1,potential%nexp_nlcc
02212              is_ok=cp_sll_val_next(list,val,error=error)
02213              CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02214              CALL val_get(val,c_val=line_att,error=error)
02215              READ(line_att,*) potential%alpha_nlcc(ipot)
02216              CALL remove_word(line_att)
02217              READ(line_att,*) potential%nct_nlcc(ipot)
02218              CALL remove_word(line_att)
02219              DO ic=1,potential%nct_nlcc(ipot)
02220                READ(line_att,*) potential%cval_nlcc(ic,ipot)
02221                !make cp2k compatible with bigdft 
02222                potential%cval_nlcc(ic,ipot)=potential%cval_nlcc(ic,ipot)/(4.0_dp*pi)
02223                CALL remove_word(line_att)
02224              END DO
02225           END DO
02226         ELSEIF(INDEX(line_att,"LSD") /= 0) THEN
02227           potential%lsdpot = .TRUE.
02228           CALL remove_word(line_att)
02229           READ(line_att,*) potential%nexp_lsd
02230           DO ipot=1,potential%nexp_lsd
02231              is_ok=cp_sll_val_next(list,val,error=error)
02232              CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02233              CALL val_get(val,c_val=line_att,error=error)
02234              READ(line_att,*) potential%alpha_lsd(ipot)
02235              CALL remove_word(line_att)
02236              READ(line_att,*) potential%nct_lsd(ipot)
02237              CALL remove_word(line_att)
02238              DO ic=1,potential%nct_lsd(ipot)
02239                READ(line_att,*) potential%cval_lsd(ic,ipot)
02240                CALL remove_word(line_att)
02241              END DO
02242           END DO
02243         ELSE
02244           EXIT
02245         END IF
02246       END DO
02247       ! Read the parameters for the non-local part of the GTH pseudopotential (ppnl)
02248       READ(line_att,*) nlmax
02249       CALL remove_word(line_att)
02250       IF (nlmax > 0) THEN
02251         ! Load the parameter for nlmax non-local projectors
02252         DO l=0,nlmax-1
02253           is_ok=cp_sll_val_next(list,val,error=error)
02254           CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02255           CALL val_get(val,c_val=line_att,error=error)
02256           READ(line_att,*) potential%rcnl(l)
02257           CALL remove_word(line_att)
02258           READ(line_att,*) potential%nl(l)
02259           CALL remove_word(line_att)
02260           DO i=1,potential%nl(l)
02261             IF (i == 1) THEN
02262               READ(line_att,*) potential%hnl(1,1,l)
02263               CALL remove_word(line_att)
02264             ELSE
02265               CPPostcondition(LEN_TRIM(line_att)==0, cp_failure_level, routineP, error, failure)
02266               is_ok=cp_sll_val_next(list,val,error=error)
02267               CPPostcondition(is_ok, cp_failure_level, routineP, error, failure)
02268               CALL val_get(val,c_val=line_att,error=error)
02269               READ(line_att,*) potential%hnl(i,i,l)
02270               CALL remove_word(line_att)
02271             END IF
02272             DO j=i+1,potential%nl(l)
02273               READ(line_att,*) potential%hnl(i,j,l)
02274               potential%hnl(j,i,l) = potential%hnl(i,j,l)
02275               CALL remove_word(line_att)
02276             END DO
02277           END DO
02278           CPPostcondition(LEN_TRIM(line_att)==0, cp_failure_level, routineP, error, failure)
02279         END DO
02280       END IF
02281     ELSE
02282       NULLIFY(parser)
02283       CALL parser_create(parser,pseudo_file,error=error)
02284 
02285       search_loop: DO
02286         CALL parser_search_string(parser,TRIM(apname),.TRUE.,found,line,error=error)
02287         IF (found) THEN
02288           CALL uppercase(symbol)
02289           CALL uppercase(apname)
02290           ! Check both the element symbol and the atomic potential name
02291           match = .FALSE.
02292           CALL uppercase(line)
02293           line2 = " "//line//" "
02294           symbol2 = " "//TRIM(symbol)//" "
02295           apname2 = " "//TRIM(apname)//" "
02296           strlen1 = LEN_TRIM(symbol2) + 1
02297           strlen2 = LEN_TRIM(apname2) + 1
02298 
02299           IF ( (INDEX(line2,symbol2(:strlen1)) > 0).AND.&
02300                   (INDEX(line2,apname2(:strlen2)) > 0) ) match = .TRUE.
02301 
02302           IF (match) THEN
02303              ! Read the electronic configuration
02304              l = 0
02305              CALL parser_get_object(parser,elec_conf(l),newline=.TRUE.,error=error)
02306              DO WHILE (parser_test_next_token(parser,error=error) == "INT")
02307                l = l + 1
02308                CALL parser_get_object(parser,elec_conf(l),error=error)
02309              END DO
02310              potential%econf(0:3) = elec_conf(0:3)
02311              potential%zion = REAL ( SUM(elec_conf), dp )
02312              ! Read r(loc) to define the exponent of the core charge
02313              CALL parser_get_object(parser,potential%rc,newline=.TRUE.,error=error)
02314              ! Read the parameters for the local part of the GTH pseudopotential (ppl)
02315              CALL parser_get_object(parser,potential%ncl,error=error)
02316              DO i=1,potential%ncl
02317                CALL parser_get_object(parser,potential%cl(i),error=error)
02318              END DO
02319              ! Extended type input
02320              DO
02321                CALL parser_get_next_line(parser,1,error=error)
02322                IF(parser_test_next_token(parser,error=error) == "INT") THEN
02323                  EXIT
02324                ELSEIF(parser_test_next_token(parser,error=error) == "STR") THEN
02325                  CALL parser_get_object(parser,line,error=error)
02326                  IF(INDEX(LINE,"LPOT") /= 0) THEN
02327                    ! local potential
02328                    potential%lpotextended = .TRUE.
02329                    CALL parser_get_object(parser,potential%nexp_lpot,error=error)
02330                    DO ipot=1,potential%nexp_lpot
02331                      CALL parser_get_object(parser,potential%alpha_lpot(ipot),newline=.TRUE.,error=error)
02332                      CALL parser_get_object(parser,potential%nct_lpot(ipot),error=error)
02333                      DO ic=1,potential%nct_lpot(ipot)
02334                        CALL parser_get_object(parser,potential%cval_lpot(ic,ipot),error=error)
02335                      END DO
02336                    END DO
02337                  ELSEIF(INDEX(LINE,"NLCC") /= 0) THEN
02338                    ! NLCC
02339                    potential%nlcc = .TRUE.
02340                    CALL parser_get_object(parser,potential%nexp_nlcc,error=error)
02341                    DO ipot=1,potential%nexp_nlcc
02342                      CALL parser_get_object(parser,potential%alpha_nlcc(ipot),newline=.TRUE.,error=error)
02343                      CALL parser_get_object(parser,potential%nct_nlcc(ipot),error=error)
02344                      DO ic=1,potential%nct_nlcc(ipot)
02345                        CALL parser_get_object(parser,potential%cval_nlcc(ic,ipot),error=error)
02346                        !make cp2k compatible with bigdft 
02347                        potential%cval_nlcc(ic,ipot)=potential%cval_nlcc(ic,ipot)/(4.0_dp*pi)
02348                      END DO
02349                    END DO
02350                  ELSEIF(INDEX(LINE,"LSD") /= 0) THEN
02351                    ! LSD potential
02352                    potential%lsdpot = .TRUE.
02353                    CALL parser_get_object(parser,potential%nexp_lsd,error=error)
02354                    DO ipot=1,potential%nexp_lsd
02355                      CALL parser_get_object(parser,potential%alpha_lsd(ipot),newline=.TRUE.,error=error)
02356                      CALL parser_get_object(parser,potential%nct_lsd(ipot),error=error)
02357                      DO ic=1,potential%nct_lsd(ipot)
02358                        CALL parser_get_object(parser,potential%cval_lsd(ic,ipot),error=error)
02359                      END DO
02360                    END DO
02361                  ELSE
02362                    CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure)
02363                  END IF
02364                ELSE
02365                  CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure)
02366                END IF
02367              END DO
02368              ! Read the parameters for the non-local part of the GTH pseudopotential (ppnl)
02369              CALL parser_get_object(parser,nlmax,error=error)
02370              IF (nlmax > 0) THEN
02371                 ! Load the parameter for n non-local projectors
02372                 DO l=0,nlmax-1
02373                    CALL parser_get_object(parser,potential%rcnl(l),newline=.TRUE.,error=error)
02374                    CALL parser_get_object(parser,potential%nl(l),error=error)
02375                    DO i=1,potential%nl(l)
02376                       IF (i == 1) THEN
02377                         CALL parser_get_object(parser,potential%hnl(i,i,l),error=error)
02378                       ELSE
02379                         CALL parser_get_object(parser,potential%hnl(i,i,l),newline=.TRUE.,error=error)
02380                       END IF
02381                       DO j=i+1,potential%nl(l)
02382                         CALL parser_get_object(parser,potential%hnl(i,j,l),error=error)
02383                         potential%hnl(j,i,l) = potential%hnl(i,j,l)
02384                       END DO
02385                    END DO
02386                 END DO
02387              END IF
02388              EXIT search_loop
02389           END IF
02390         ELSE
02391           ! Stop program, if the end of file is reached
02392           CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure)
02393         END IF
02394 
02395       END DO search_loop
02396 
02397       CALL parser_release(parser,error=error)
02398     END IF
02399 
02400   END SUBROUTINE read_gth_potential
02401 ! *****************************************************************************
02402 
02403 END MODULE atom_types