CP2K 2.4 (Revision 12889)

basis_set_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 ! *****************************************************************************
00011 MODULE basis_set_types
00012 
00013   USE ai_coulomb,                      ONLY: coulomb2
00014   USE bibliography,                    ONLY: VandeVondele2007,&
00015                                              cite_reference
00016   USE cp_linked_list_val,              ONLY: cp_sll_val_next,&
00017                                              cp_sll_val_type
00018   USE cp_para_types,                   ONLY: cp_para_env_type
00019   USE cp_parser_methods,               ONLY: parser_get_object,&
00020                                              parser_search_string
00021   USE cp_parser_types,                 ONLY: cp_parser_type,&
00022                                              parser_create,&
00023                                              parser_release
00024   USE input_section_types,             ONLY: section_vals_get,&
00025                                              section_vals_list_get,&
00026                                              section_vals_type,&
00027                                              section_vals_val_get,&
00028                                              section_vals_val_set
00029   USE input_val_types,                 ONLY: val_get,&
00030                                              val_type
00031   USE kinds,                           ONLY: default_path_length,&
00032                                              default_string_length,&
00033                                              dp
00034   USE mathconstants,                   ONLY: dfac,&
00035                                              pi
00036   USE memory_utilities,                ONLY: reallocate
00037   USE orbital_pointers,                ONLY: coset,&
00038                                              indco,&
00039                                              init_orbital_pointers,&
00040                                              nco,&
00041                                              ncoset,&
00042                                              nso,&
00043                                              nsoset
00044   USE orbital_symbols,                 ONLY: cgf_symbol,&
00045                                              sgf_symbol
00046   USE orbital_transformation_matrices, ONLY: orbtramat
00047   USE sto_ng,                          ONLY: get_sto_ng
00048   USE string_utilities,                ONLY: integer_to_string,&
00049                                              remove_word,&
00050                                              uppercase
00051   USE termination,                     ONLY: stop_program
00052   USE timings,                         ONLY: timeset,&
00053                                              timestop
00054 #include "cp_common_uses.h"
00055 
00056   IMPLICIT NONE
00057 
00058   PRIVATE
00059 
00060   ! Global parameters (only in this module)
00061 
00062   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'basis_set_types'
00063 
00064 ! *****************************************************************************
00065   ! Define the Gaussian-type orbital basis set type
00066 
00067   TYPE gto_basis_set_type
00068      !MK PRIVATE
00069      CHARACTER(LEN=default_string_length)       :: name
00070      REAL(KIND = dp)                            :: kind_radius
00071      REAL(KIND = dp)                            :: short_kind_radius
00072      INTEGER                                    :: norm_type
00073      INTEGER                                    :: ncgf,nset,nsgf
00074      CHARACTER(LEN=12), DIMENSION(:), POINTER   :: cgf_symbol
00075      CHARACTER(LEN=6), DIMENSION(:), POINTER    :: sgf_symbol
00076      REAL(KIND = dp), DIMENSION(:), POINTER     :: norm_cgf,set_radius
00077      INTEGER, DIMENSION(:), POINTER             :: lmax,lmin,lx,ly,lz,m,ncgf_set,
00078           npgf,nsgf_set,nshell
00079      REAL(KIND = dp), DIMENSION(:,:), POINTER   :: cphi,pgf_radius,sphi,zet
00080      INTEGER, DIMENSION(:,:), POINTER           :: first_cgf,first_sgf,l,
00081           last_cgf,last_sgf,n
00082      REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: gcc
00083   END TYPE gto_basis_set_type
00084 
00085   TYPE gto_basis_set_p_type
00086      TYPE(gto_basis_set_type), POINTER :: gto_basis_set
00087   END TYPE gto_basis_set_p_type
00088 
00089 ! *****************************************************************************
00090   ! Define the Slater-type orbital basis set type
00091 
00092   TYPE sto_basis_set_type
00093      !MK PRIVATE
00094      CHARACTER(LEN=default_string_length)       :: name
00095      INTEGER                                    :: nshell
00096      CHARACTER(LEN=6), DIMENSION(:), POINTER    :: symbol
00097      INTEGER, DIMENSION(:), POINTER             :: nq,lq
00098      REAL(KIND = dp), DIMENSION(:), POINTER     :: zet
00099   END TYPE sto_basis_set_type
00100 
00101   TYPE sto_basis_set_p_type
00102      TYPE(sto_basis_set_type), POINTER :: sto_basis_set
00103   END TYPE sto_basis_set_p_type
00104 
00105 ! *****************************************************************************
00106   ! Define the Geminal basis set type
00107 
00108   TYPE geminal_basis_set_type
00109      CHARACTER(LEN=default_string_length)       :: name
00110      CHARACTER(LEN=12), DIMENSION(:), POINTER   :: cgf_symbol
00111      CHARACTER(LEN=2)                           :: type_restriction
00112      REAL(KIND = dp)                            :: kind_radius
00113      REAL(KIND = dp), DIMENSION(:), POINTER     :: set_radius
00114      REAL(KIND = dp), DIMENSION(:,:), POINTER   :: pgf_radius
00115      INTEGER                                    :: nset
00116      INTEGER                                    :: ngeminals
00117      INTEGER, DIMENSION(:), POINTER             :: lmax,lmin,ls
00118      INTEGER, DIMENSION(:), POINTER             :: npgf
00119      INTEGER, DIMENSION(:), POINTER             :: ngem_set,nshell
00120      INTEGER, DIMENSION(:,:), POINTER           :: l
00121      REAL(KIND = dp), DIMENSION(:,:,:,:), 
00122                                         POINTER :: zet, zeth
00123      INTEGER, DIMENSION(:,:), POINTER           :: first_cgf,last_cgf
00124      REAL(KIND = dp), DIMENSION(:,:,:), POINTER :: gcc
00125   END TYPE geminal_basis_set_type
00126 
00127   TYPE geminal_basis_set_p_type
00128      TYPE(geminal_basis_set_type), POINTER :: geminal_basis_set
00129   END TYPE geminal_basis_set_p_type
00130 
00131 ! *****************************************************************************
00132 
00133   ! Public subroutines
00134   PUBLIC :: allocate_gto_basis_set,&
00135             deallocate_gto_basis_set,&
00136             get_gto_basis_set_info,&
00137             get_gto_basis_set,&
00138             init_aux_basis_set,&
00139             init_cphi_and_sphi,&
00140             init_orb_basis_set,&
00141             read_gto_basis_set,&
00142             copy_gto_basis_set,&
00143             set_gto_basis_set,&
00144             write_gto_basis_set,&
00145             write_orb_basis_set
00146 
00147   PUBLIC :: allocate_sto_basis_set,&
00148             create_gto_from_sto_basis,&
00149             deallocate_sto_basis_set,&
00150             get_sto_basis_set,&
00151             set_sto_basis_set,&
00152             srules
00153 
00154   PUBLIC :: allocate_geminal_basis_set,&
00155             deallocate_geminal_basis_set,&
00156             get_geminal_basis_set,&
00157             read_geminal_basis_set,&
00158             set_geminal_basis_set,&
00159             write_geminal_basis_set
00160 
00161   ! Public data types
00162   PUBLIC :: gto_basis_set_p_type,&
00163             gto_basis_set_type,&
00164             sto_basis_set_p_type,&
00165             sto_basis_set_type,&
00166             geminal_basis_set_p_type,&
00167             geminal_basis_set_type
00168 
00169 CONTAINS
00170 
00171 ! *****************************************************************************
00172   SUBROUTINE allocate_gto_basis_set(gto_basis_set, error)
00173 
00174     ! Allocate a Gaussian-type orbital (GTO) basis set data set.
00175 
00176     ! - Creation (26.10.2000,MK)
00177 
00178     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00179     TYPE(cp_error_type), INTENT(inout)       :: error
00180 
00181     CHARACTER(len=*), PARAMETER :: routineN = 'allocate_gto_basis_set', 
00182       routineP = moduleN//':'//routineN
00183 
00184     INTEGER                                  :: stat
00185     LOGICAL                                  :: failure
00186 
00187     failure = .FALSE.
00188     CALL deallocate_gto_basis_set(gto_basis_set,error)
00189 
00190     ALLOCATE (gto_basis_set,STAT=stat)
00191     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00192 
00193     NULLIFY (gto_basis_set%cgf_symbol)
00194     NULLIFY (gto_basis_set%first_cgf)
00195     NULLIFY (gto_basis_set%first_sgf)
00196     NULLIFY (gto_basis_set%gcc)
00197     NULLIFY (gto_basis_set%l)
00198     NULLIFY (gto_basis_set%last_cgf)
00199     NULLIFY (gto_basis_set%last_sgf)
00200     NULLIFY (gto_basis_set%lmax)
00201     NULLIFY (gto_basis_set%lmin)
00202     NULLIFY (gto_basis_set%lx)
00203     NULLIFY (gto_basis_set%ly)
00204     NULLIFY (gto_basis_set%lz)
00205     NULLIFY (gto_basis_set%m)
00206     NULLIFY (gto_basis_set%n)
00207     NULLIFY (gto_basis_set%ncgf_set)
00208     NULLIFY (gto_basis_set%norm_cgf)
00209     NULLIFY (gto_basis_set%npgf)
00210     NULLIFY (gto_basis_set%nsgf_set)
00211     NULLIFY (gto_basis_set%nshell)
00212     NULLIFY (gto_basis_set%pgf_radius)
00213     NULLIFY (gto_basis_set%cphi)
00214     NULLIFY (gto_basis_set%sphi)
00215     NULLIFY (gto_basis_set%set_radius)
00216     NULLIFY (gto_basis_set%sgf_symbol)
00217     NULLIFY (gto_basis_set%zet)
00218     gto_basis_set%short_kind_radius = 0.0_dp
00219 
00220   END SUBROUTINE allocate_gto_basis_set
00221 
00222 ! *****************************************************************************
00223   SUBROUTINE deallocate_gto_basis_set(gto_basis_set, error)
00224 
00225     ! Deallocate a Gaussian-type orbital (GTO) basis set data set.
00226 
00227     ! - Creation (03.11.2000,MK)
00228 
00229     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00230     TYPE(cp_error_type), INTENT(inout)       :: error
00231 
00232     CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_gto_basis_set', 
00233       routineP = moduleN//':'//routineN
00234 
00235     INTEGER                                  :: stat
00236     LOGICAL                                  :: failure
00237 
00238     failure = .FALSE.
00239     IF (ASSOCIATED(gto_basis_set))  THEN
00240        IF (ASSOCIATED(gto_basis_set%cgf_symbol)) THEN
00241           DEALLOCATE (gto_basis_set%cgf_symbol,STAT=stat)
00242           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00243        ENDIF
00244        IF (ASSOCIATED(gto_basis_set%sgf_symbol)) THEN
00245           DEALLOCATE (gto_basis_set%sgf_symbol,STAT=stat)
00246           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00247        ENDIF
00248        DEALLOCATE (gto_basis_set%norm_cgf,STAT=stat)
00249        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00250        DEALLOCATE (gto_basis_set%set_radius,STAT=stat)
00251        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00252        DEALLOCATE (gto_basis_set%lmax,STAT=stat)
00253        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00254        DEALLOCATE (gto_basis_set%lmin,STAT=stat)
00255        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00256        DEALLOCATE (gto_basis_set%lx,STAT=stat)
00257        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00258        DEALLOCATE (gto_basis_set%ly,STAT=stat)
00259        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00260        DEALLOCATE (gto_basis_set%lz,STAT=stat)
00261        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00262        DEALLOCATE (gto_basis_set%m,STAT=stat)
00263        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00264        DEALLOCATE (gto_basis_set%ncgf_set,STAT=stat)
00265        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00266        DEALLOCATE (gto_basis_set%npgf,STAT=stat)
00267        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00268        DEALLOCATE (gto_basis_set%nsgf_set,STAT=stat)
00269        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00270        DEALLOCATE (gto_basis_set%nshell,STAT=stat)
00271        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00272        DEALLOCATE (gto_basis_set%cphi,STAT=stat)
00273        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00274        DEALLOCATE (gto_basis_set%pgf_radius,STAT=stat)
00275        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00276        DEALLOCATE (gto_basis_set%sphi,STAT=stat)
00277        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00278        DEALLOCATE (gto_basis_set%zet,STAT=stat)
00279        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00280        DEALLOCATE (gto_basis_set%first_cgf,STAT=stat)
00281        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00282        DEALLOCATE (gto_basis_set%first_sgf,STAT=stat)
00283        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00284        DEALLOCATE (gto_basis_set%l,STAT=stat)
00285        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00286        DEALLOCATE (gto_basis_set%last_cgf,STAT=stat)
00287        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00288        DEALLOCATE (gto_basis_set%last_sgf,STAT=stat)
00289        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00290        DEALLOCATE (gto_basis_set%n,STAT=stat)
00291        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00292        DEALLOCATE (gto_basis_set%gcc,STAT=stat)
00293        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00294        DEALLOCATE (gto_basis_set,STAT=stat)
00295        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00296     END IF
00297   END SUBROUTINE deallocate_gto_basis_set
00298 
00299 ! *****************************************************************************
00300   SUBROUTINE copy_gto_basis_set(basis_set_in, basis_set_out, error)
00301 
00302     ! Copy a Gaussian-type orbital (GTO) basis set data set.
00303 
00304     TYPE(gto_basis_set_type), POINTER        :: basis_set_in, basis_set_out
00305     TYPE(cp_error_type), INTENT(inout)       :: error
00306 
00307     CHARACTER(len=*), PARAMETER :: routineN = 'copy_gto_basis_set', 
00308       routineP = moduleN//':'//routineN
00309 
00310     INTEGER                                  :: maxco, maxpgf, maxshell, 
00311                                                 ncgf, nset, nsgf, stat
00312     LOGICAL                                  :: failure
00313     TYPE(gto_basis_set_type), POINTER        :: bin, bout
00314 
00315     failure = .FALSE.
00316 
00317     CALL deallocate_gto_basis_set(basis_set_out,error)
00318     ALLOCATE (basis_set_out,STAT=stat)
00319     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00320     CPPrecondition(ASSOCIATED(basis_set_in),cp_failure_level,routineP,error,failure)
00321 
00322     bin => basis_set_in
00323     bout => basis_set_out
00324 
00325     bout%name = bin%name
00326     bout%kind_radius = bin%kind_radius
00327     bout%norm_type = bin%norm_type
00328     bout%nset = bin%nset
00329     bout%ncgf = bin%ncgf
00330     bout%nsgf = bin%nsgf
00331     nset = bin%nset
00332     ncgf = bin%ncgf
00333     nsgf = bin%nsgf
00334     ALLOCATE (bout%cgf_symbol(ncgf),STAT=stat)
00335     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00336     ALLOCATE (bout%sgf_symbol(nsgf),STAT=stat)
00337     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00338     bout%cgf_symbol = bin%cgf_symbol
00339     bout%sgf_symbol = bin%sgf_symbol
00340     ALLOCATE (bout%norm_cgf(ncgf),STAT=stat)
00341     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00342     bout%norm_cgf = bin%norm_cgf
00343     ALLOCATE (bout%set_radius(nset),STAT=stat)
00344     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00345     bout%set_radius = bin%set_radius
00346     ALLOCATE (bout%lmax(nset),bout%lmin(nset),bout%npgf(nset),bout%nshell(nset),STAT=stat)
00347     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00348     bout%lmax = bin%lmax
00349     bout%lmin = bin%lmin
00350     bout%npgf = bin%npgf
00351     bout%nshell = bin%nshell
00352     ALLOCATE (bout%lx(ncgf),bout%ly(ncgf),bout%lz(ncgf),bout%m(nsgf),STAT=stat)
00353     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00354     bout%lx = bin%lx
00355     bout%ly = bin%ly
00356     bout%lz = bin%lz
00357     bout%m = bin%m
00358     ALLOCATE (bout%ncgf_set(nset),bout%nsgf_set(nset),STAT=stat)
00359     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00360     bout%ncgf_set = bin%ncgf_set
00361     bout%nsgf_set = bin%nsgf_set
00362     maxco = SIZE(bin%cphi,1)
00363     ALLOCATE (bout%cphi(maxco,ncgf),bout%sphi(maxco,nsgf),STAT=stat)
00364     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00365     bout%cphi = bin%cphi
00366     bout%sphi = bin%sphi
00367     maxpgf = MAXVAL(bin%npgf)
00368     ALLOCATE (bout%pgf_radius(maxpgf,nset),bout%zet(maxpgf,nset),STAT=stat)
00369     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00370     bout%pgf_radius = bin%pgf_radius
00371     bout%zet = bin%zet
00372     maxshell = MAXVAL(bin%nshell)
00373     ALLOCATE (bout%first_cgf(maxshell,nset),bout%first_sgf(maxshell,nset),STAT=stat)
00374     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00375     ALLOCATE (bout%last_cgf(maxshell,nset),bout%last_sgf(maxshell,nset),STAT=stat)
00376     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00377     bout%first_cgf = bin%first_cgf
00378     bout%first_sgf = bin%first_sgf
00379     bout%last_cgf = bin%last_cgf
00380     bout%last_sgf = bin%last_sgf
00381     ALLOCATE (bout%n(maxshell,nset),bout%l(maxshell,nset),STAT=stat)
00382     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00383     bout%n = bin%n
00384     bout%l = bin%l
00385     ALLOCATE (bout%gcc(maxpgf,maxshell,nset),STAT=stat)
00386     CPPostcondition(stat==0, cp_failure_level, routineP, error, failure)
00387     bout%gcc = bin%gcc
00388 
00389   END SUBROUTINE copy_gto_basis_set
00390 
00391 ! *****************************************************************************
00392   SUBROUTINE get_gto_basis_set_info(gto_basis,maxco,maxl,maxsgf,maxder)
00393 
00394     TYPE(gto_basis_set_p_type), 
00395       DIMENSION(:), POINTER                  :: gto_basis
00396     INTEGER, INTENT(OUT), OPTIONAL           :: maxco, maxl, maxsgf
00397     INTEGER, INTENT(IN), OPTIONAL            :: maxder
00398 
00399     CHARACTER(len=*), PARAMETER :: routineN = 'get_gto_basis_set_info', 
00400       routineP = moduleN//':'//routineN
00401 
00402     INTEGER                                  :: ibas, nbas, val
00403     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00404 
00405 ! -------------------------------------------------------------------------
00406 
00407     IF (ASSOCIATED(gto_basis)) THEN
00408        nbas = SIZE(gto_basis)
00409        IF (PRESENT(maxco)) THEN
00410           maxco = 0
00411           DO ibas=1,nbas
00412              gto_basis_set => gto_basis(ibas)%gto_basis_set
00413              IF (PRESENT(maxder)) THEN
00414                 CALL get_gto_basis_set(gto_basis_set,maxco=val,maxder=maxder)
00415              ELSE
00416                 CALL get_gto_basis_set(gto_basis_set,maxco=val)
00417              END IF
00418              maxco = MAX(val,maxco)
00419           END DO
00420        END IF
00421        IF (PRESENT(maxl)) THEN
00422           maxl = 0
00423           DO ibas=1,nbas
00424              gto_basis_set => gto_basis(ibas)%gto_basis_set
00425              CALL get_gto_basis_set(gto_basis_set,maxl=val)
00426              maxl = MAX(val,maxl)
00427           END DO
00428        END IF
00429        IF (PRESENT(maxsgf)) THEN
00430           maxsgf = 0
00431           DO ibas=1,nbas
00432              gto_basis_set => gto_basis(ibas)%gto_basis_set
00433              CALL get_gto_basis_set(gto_basis_set,nsgf=val)
00434              maxsgf = MAX(val,maxsgf)
00435           END DO
00436        END IF
00437     ELSE
00438        CALL stop_program(routineN,moduleN,__LINE__,&
00439             "The pointer gto_basis is not associated")
00440     END IF
00441 
00442   END SUBROUTINE get_gto_basis_set_info
00443 ! *****************************************************************************
00444   SUBROUTINE get_gto_basis_set(gto_basis_set,name,norm_type,kind_radius,ncgf,&
00445        nset,nsgf,cgf_symbol,sgf_symbol,norm_cgf,set_radius,lmax,lmin,lx,ly,lz,&
00446        m,ncgf_set,npgf,nsgf_set,nshell,cphi,pgf_radius,sphi,zet,first_cgf,first_sgf,l,&
00447        last_cgf,last_sgf,n,gcc,maxco,maxl,maxpgf,maxsgf_set,maxshell,maxso,nco_sum,&
00448        npgf_sum,nshell_sum,maxder,short_kind_radius)
00449 
00450     ! Get informations about a Gaussian-type orbital (GTO) basis set.
00451 
00452     ! - Creation (10.01.2002,MK)
00453 
00454     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00455     CHARACTER(LEN=default_string_length), 
00456       INTENT(OUT), OPTIONAL                  :: name
00457     INTEGER, INTENT(OUT), OPTIONAL           :: norm_type
00458     REAL(KIND=dp), INTENT(OUT), OPTIONAL     :: kind_radius
00459     INTEGER, INTENT(OUT), OPTIONAL           :: ncgf, nset, nsgf
00460     CHARACTER(LEN=12), DIMENSION(:), 
00461       OPTIONAL, POINTER                      :: cgf_symbol
00462     CHARACTER(LEN=6), DIMENSION(:), 
00463       OPTIONAL, POINTER                      :: sgf_symbol
00464     REAL(KIND=dp), DIMENSION(:), OPTIONAL, 
00465       POINTER                                :: norm_cgf, set_radius
00466     INTEGER, DIMENSION(:), OPTIONAL, POINTER :: lmax, lmin, lx, ly, lz, m, 
00467                                                 ncgf_set, npgf, nsgf_set, 
00468                                                 nshell
00469     REAL(KIND=dp), DIMENSION(:, :), 
00470       OPTIONAL, POINTER                      :: cphi, pgf_radius, sphi, zet
00471     INTEGER, DIMENSION(:, :), OPTIONAL, 
00472       POINTER                                :: first_cgf, first_sgf, l, 
00473                                                 last_cgf, last_sgf, n
00474     REAL(KIND=dp), DIMENSION(:, :, :), 
00475       OPTIONAL, POINTER                      :: gcc
00476     INTEGER, INTENT(OUT), OPTIONAL           :: maxco, maxl, maxpgf, 
00477                                                 maxsgf_set, maxshell, maxso, 
00478                                                 nco_sum, npgf_sum, nshell_sum
00479     INTEGER, INTENT(IN), OPTIONAL            :: maxder
00480     REAL(KIND=dp), INTENT(OUT), OPTIONAL     :: short_kind_radius
00481 
00482     CHARACTER(len=*), PARAMETER :: routineN = 'get_gto_basis_set', 
00483       routineP = moduleN//':'//routineN
00484 
00485     INTEGER                                  :: iset, nder
00486 
00487 ! -------------------------------------------------------------------------
00488 
00489     IF (ASSOCIATED(gto_basis_set)) THEN
00490        IF (PRESENT(name)) name = gto_basis_set%name
00491        IF (PRESENT(norm_type)) norm_type = gto_basis_set%norm_type
00492        IF (PRESENT(kind_radius)) kind_radius = gto_basis_set%kind_radius
00493        IF (PRESENT(short_kind_radius)) short_kind_radius = gto_basis_set%short_kind_radius
00494        IF (PRESENT(ncgf)) ncgf = gto_basis_set%ncgf
00495        IF (PRESENT(nset)) nset = gto_basis_set%nset
00496        IF (PRESENT(nsgf)) nsgf = gto_basis_set%nsgf
00497        IF (PRESENT(cgf_symbol)) cgf_symbol => gto_basis_set%cgf_symbol
00498        IF (PRESENT(sgf_symbol)) sgf_symbol => gto_basis_set%sgf_symbol
00499        IF (PRESENT(norm_cgf)) norm_cgf => gto_basis_set%norm_cgf
00500        IF (PRESENT(set_radius)) set_radius => gto_basis_set%set_radius
00501        IF (PRESENT(lmax)) lmax => gto_basis_set%lmax
00502        IF (PRESENT(lmin)) lmin => gto_basis_set%lmin
00503        IF (PRESENT(lx)) lx => gto_basis_set%lx
00504        IF (PRESENT(ly)) ly => gto_basis_set%ly
00505        IF (PRESENT(lz)) lz => gto_basis_set%lz
00506        IF (PRESENT(m)) m => gto_basis_set%m
00507        IF (PRESENT(ncgf_set)) ncgf_set => gto_basis_set%ncgf_set
00508        IF (PRESENT(npgf)) npgf => gto_basis_set%npgf
00509        IF (PRESENT(nsgf_set)) nsgf_set => gto_basis_set%nsgf_set
00510        IF (PRESENT(nshell)) nshell => gto_basis_set%nshell
00511        IF (PRESENT(cphi)) cphi => gto_basis_set%cphi
00512        IF (PRESENT(pgf_radius)) pgf_radius => gto_basis_set%pgf_radius
00513        IF (PRESENT(sphi)) sphi => gto_basis_set%sphi
00514        IF (PRESENT(zet)) zet => gto_basis_set%zet
00515        IF (PRESENT(first_cgf)) first_cgf => gto_basis_set%first_cgf
00516        IF (PRESENT(first_sgf)) first_sgf => gto_basis_set%first_sgf
00517        IF (PRESENT(l)) l => gto_basis_set%l
00518        IF (PRESENT(last_cgf)) last_cgf => gto_basis_set%last_cgf
00519        IF (PRESENT(last_sgf)) last_sgf => gto_basis_set%last_sgf
00520        IF (PRESENT(n)) n => gto_basis_set%n
00521        IF (PRESENT(gcc)) gcc => gto_basis_set%gcc
00522        IF (PRESENT(maxco)) THEN
00523           maxco = 0
00524           IF (PRESENT(maxder)) THEN
00525              nder = maxder
00526           ELSE
00527              nder = 0
00528           END IF
00529           DO iset=1,gto_basis_set%nset
00530              maxco = MAX(maxco,gto_basis_set%npgf(iset)*&
00531                   ncoset(gto_basis_set%lmax(iset)+nder))
00532           END DO
00533        END IF
00534        IF (PRESENT(maxl)) THEN
00535           maxl = -1
00536           DO iset=1,gto_basis_set%nset
00537              maxl = MAX(maxl,gto_basis_set%lmax(iset))
00538           END DO
00539        END IF
00540        IF (PRESENT(maxpgf)) THEN
00541           maxpgf = 0
00542           DO iset=1,gto_basis_set%nset
00543              maxpgf = MAX(maxpgf,gto_basis_set%npgf(iset))
00544           END DO
00545        END IF
00546        IF (PRESENT(maxsgf_set)) THEN
00547           maxsgf_set = 0
00548           DO iset=1,gto_basis_set%nset
00549              maxsgf_set = MAX(maxsgf_set,gto_basis_set%nsgf_set(iset))
00550           END DO
00551        END IF
00552        IF (PRESENT(maxshell)) THEN ! MAXVAL on structure component avoided
00553           maxshell = 0
00554           DO iset=1,gto_basis_set%nset
00555              maxshell = MAX(maxshell,gto_basis_set%nshell(iset))
00556           END DO
00557        END IF
00558        IF (PRESENT(maxso)) THEN
00559           maxso = 0
00560           DO iset=1,gto_basis_set%nset
00561              maxso = MAX(maxso,gto_basis_set%npgf(iset)*&
00562                   nsoset(gto_basis_set%lmax(iset)))
00563           END DO
00564        END IF
00565 
00566        IF (PRESENT(nco_sum)) THEN
00567           nco_sum = 0
00568           DO iset=1,gto_basis_set%nset
00569              nco_sum = nco_sum + gto_basis_set%npgf(iset)*&
00570                   ncoset(gto_basis_set%lmax(iset))
00571           END DO
00572        END IF
00573        IF (PRESENT(npgf_sum)) npgf_sum = SUM(gto_basis_set%npgf)
00574        IF (PRESENT(nshell_sum)) nshell_sum = SUM(gto_basis_set%nshell)
00575     ELSE
00576        CALL stop_program(routineN,moduleN,__LINE__,&
00577             "The pointer gto_basis_set is not associated")
00578     END IF
00579 
00580   END SUBROUTINE get_gto_basis_set
00581 
00582 ! *****************************************************************************
00583   SUBROUTINE init_aux_basis_set(gto_basis_set,error)
00584 
00585     ! Initialise a Gaussian-type orbital (GTO) basis set data set.
00586 
00587     ! - Creation (06.12.2000,MK)
00588 
00589     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00590     TYPE(cp_error_type), INTENT(inout)       :: error
00591 
00592     CHARACTER(len=*), PARAMETER :: routineN = 'init_aux_basis_set', 
00593       routineP = moduleN//':'//routineN
00594 
00595     INTEGER                                  :: handle
00596 
00597 ! -------------------------------------------------------------------------
00598 
00599     IF (.NOT.ASSOCIATED(gto_basis_set)) RETURN
00600 
00601     CALL timeset(routineN,handle)
00602 
00603     SELECT CASE (gto_basis_set%norm_type)
00604     CASE ( 0 )
00605        ! No normalisation requested
00606     CASE ( 1 )
00607        CALL init_norm_cgf_aux_2(gto_basis_set,error)
00608     CASE ( 2 )
00609        ! WARNING this was never tested
00610        CALL init_norm_cgf_aux(gto_basis_set,error)
00611     CASE DEFAULT
00612        CALL stop_program(routineN,moduleN,__LINE__,&
00613             "Normalization method not specified")
00614     END SELECT
00615 
00616     ! Initialise the transformation matrices "pgf" -> "cgf"
00617     CALL init_cphi_and_sphi(gto_basis_set,error)
00618   
00619     CALL timestop(handle)
00620 
00621   END SUBROUTINE init_aux_basis_set
00622 
00623 ! *****************************************************************************
00624   SUBROUTINE init_cphi_and_sphi(gto_basis_set,error)
00625 
00626     ! Initialise the matrices for the transformation of primitive Cartesian
00627     ! Gaussian-type functions to contracted Cartesian (cphi) and spherical
00628     ! (sphi) Gaussian-type functions.
00629 
00630     ! - Creation (20.09.2000,MK)
00631 
00632     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00633     TYPE(cp_error_type), INTENT(inout)       :: error
00634 
00635     CHARACTER(len=*), PARAMETER :: routineN = 'init_cphi_and_sphi', 
00636       routineP = moduleN//':'//routineN
00637 
00638     INTEGER                                  :: first_cgf, first_sgf, handle, 
00639                                                 icgf, ico, ipgf, iset, 
00640                                                 ishell, l, n, ncgf, nsgf
00641 
00642 ! -------------------------------------------------------------------------
00643 ! Build the Cartesian transformation matrix "cphi"
00644 
00645     CALL timeset(routineN,handle)
00646 
00647     DO iset=1,gto_basis_set%nset
00648        n = ncoset(gto_basis_set%lmax(iset))
00649        DO ishell=1,gto_basis_set%nshell(iset)
00650           DO icgf=gto_basis_set%first_cgf(ishell,iset),&
00651                gto_basis_set%last_cgf(ishell,iset)
00652              ico = coset(gto_basis_set%lx(icgf),&
00653                   gto_basis_set%ly(icgf),&
00654                   gto_basis_set%lz(icgf))
00655              DO ipgf=1,gto_basis_set%npgf(iset)
00656                 gto_basis_set%cphi(ico,icgf) = gto_basis_set%norm_cgf(icgf)*&
00657                      gto_basis_set%gcc(ipgf,ishell,iset)
00658                 ico = ico + n
00659              END DO
00660           END DO
00661        END DO
00662     END DO
00663 
00664     ! Build the spherical transformation matrix "sphi"
00665 
00666     n = SIZE(gto_basis_set%cphi,1)
00667 
00668     IF(n.GT.0) THEN
00669        DO iset=1,gto_basis_set%nset
00670           DO ishell=1,gto_basis_set%nshell(iset)
00671              l = gto_basis_set%l(ishell,iset)
00672              first_cgf = gto_basis_set%first_cgf(ishell,iset)
00673              first_sgf = gto_basis_set%first_sgf(ishell,iset)
00674              ncgf = nco(l)
00675              nsgf = nso(l)
00676              CALL dgemm("N","T",n,nsgf,ncgf,&
00677                   1.0_dp,gto_basis_set%cphi(1,first_cgf),n,&
00678                   orbtramat(l)%c2s(1,1),nsgf,&
00679                   0.0_dp,gto_basis_set%sphi(1,first_sgf),n)
00680           END DO
00681        END DO
00682     ENDIF
00683 
00684     CALL timestop(handle)
00685 
00686   END SUBROUTINE init_cphi_and_sphi
00687 
00688 ! *****************************************************************************
00689   SUBROUTINE init_norm_cgf_aux(gto_basis_set,error)
00690 
00691     ! Initialise the normalization factors of the contracted Cartesian Gaussian
00692     ! functions, if the Gaussian functions represent charge distributions.
00693 
00694     ! - Creation (07.12.2000,MK)
00695 
00696     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00697     TYPE(cp_error_type), INTENT(inout)       :: error
00698 
00699     CHARACTER(len=*), PARAMETER :: routineN = 'init_norm_cgf_aux', 
00700       routineP = moduleN//':'//routineN
00701 
00702     INTEGER                                  :: icgf, ico, ipgf, iset, 
00703                                                 ishell, jco, jpgf, ll, lmax, 
00704                                                 lmin, lx, ly, lz, n, npgfa, 
00705                                                 stat
00706     LOGICAL                                  :: failure
00707     REAL(KIND=dp)                            :: fnorm, gcca, gccb
00708     REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ff
00709     REAL(KIND=dp), ALLOCATABLE, 
00710       DIMENSION(:, :)                        :: gaa
00711     REAL(KIND=dp), ALLOCATABLE, 
00712       DIMENSION(:, :, :)                     :: vv
00713     REAL(KIND=dp), DIMENSION(:), POINTER     :: rpgfa, zeta
00714 
00715 ! -------------------------------------------------------------------------
00716 
00717     failure = .FALSE.
00718     n = 0
00719     ll = 0
00720     DO iset=1,gto_basis_set%nset
00721        n = MAX(n,gto_basis_set%npgf(iset)*ncoset(gto_basis_set%lmax(iset)))
00722        ll = MAX(ll,gto_basis_set%lmax(iset))
00723     END DO
00724 
00725     ALLOCATE (gaa(n,n),STAT=stat)
00726     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00727     ALLOCATE (vv(ncoset(ll),ncoset(ll),ll+ll+1),STAT=stat)
00728     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00729     ALLOCATE (ff(0:ll+ll),STAT=stat)
00730     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00731 
00732     DO iset=1,gto_basis_set%nset
00733        lmax = gto_basis_set%lmax(iset)
00734        lmin = gto_basis_set%lmin(iset)
00735        n = ncoset(lmax)
00736        npgfa = gto_basis_set%npgf(iset)
00737        rpgfa => gto_basis_set%pgf_radius(1:npgfa,iset)
00738        zeta => gto_basis_set%zet(1:npgfa,iset)
00739        CALL coulomb2(lmax,npgfa,zeta,rpgfa,lmin,&
00740             lmax,npgfa,zeta,rpgfa,lmin,&
00741             (/0.0_dp,0.0_dp,0.0_dp/),0.0_dp,gaa,vv,ff(0:))
00742        DO ishell=1,gto_basis_set%nshell(iset)
00743           DO icgf=gto_basis_set%first_cgf(ishell,iset),&
00744                gto_basis_set%last_cgf(ishell,iset)
00745              lx = gto_basis_set%lx(icgf)
00746              ly = gto_basis_set%ly(icgf)
00747              lz = gto_basis_set%lz(icgf)
00748              ico = coset(lx,ly,lz)
00749              fnorm = 0.0_dp
00750              DO ipgf=1,npgfa
00751                 gcca = gto_basis_set%gcc(ipgf,ishell,iset)
00752                 jco = coset(lx,ly,lz)
00753                 DO jpgf=1,npgfa
00754                    gccb = gto_basis_set%gcc(jpgf,ishell,iset)
00755                    fnorm = fnorm + gcca*gccb*gaa(ico,jco)
00756                    jco = jco + n
00757                 END DO
00758                 ico = ico + n
00759              END DO
00760              gto_basis_set%norm_cgf(icgf) = 1.0_dp/SQRT(fnorm)
00761           END DO
00762        END DO
00763     END DO
00764 
00765     DEALLOCATE (vv,ff, STAT=stat)
00766     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00767     DEALLOCATE (gaa,STAT=stat)
00768     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
00769 
00770   END SUBROUTINE init_norm_cgf_aux
00771 
00772 ! *****************************************************************************
00773   SUBROUTINE init_norm_cgf_aux_2(gto_basis_set,error)
00774 
00775     ! Initialise the normalization factors of the auxiliary Cartesian Gaussian
00776     ! functions (Kim-Gordon polarization basis) Norm = 1.
00777 
00778     ! - Creation (07.12.2000,GT)
00779 
00780     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00781     TYPE(cp_error_type), INTENT(inout)       :: error
00782 
00783     INTEGER                                  :: icgf, iset, ishell
00784 
00785 ! -------------------------------------------------------------------------
00786 
00787     DO iset=1,gto_basis_set%nset
00788        DO ishell=1,gto_basis_set%nshell(iset)
00789           DO icgf=gto_basis_set%first_cgf(ishell,iset),&
00790                gto_basis_set%last_cgf(ishell,iset)
00791              gto_basis_set%norm_cgf(icgf) = 1.0_dp
00792           END DO
00793        END DO
00794     END DO
00795 
00796   END SUBROUTINE init_norm_cgf_aux_2
00797 
00798 ! *****************************************************************************
00799   SUBROUTINE init_norm_cgf_orb(gto_basis_set,error)
00800 
00801     ! Initialise the normalization factors of the contracted Cartesian Gaussian
00802     ! functions.
00803 
00804     ! - Creation (14.04.2000,MK)
00805 
00806     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00807     TYPE(cp_error_type), INTENT(inout)       :: error
00808 
00809     INTEGER                                  :: icgf, ipgf, iset, ishell, 
00810                                                 jpgf, l, lx, ly, lz
00811     REAL(KIND=dp)                            :: expzet, fnorm, gcca, gccb, 
00812                                                 prefac, zeta, zetb
00813 
00814 ! -------------------------------------------------------------------------
00815 
00816     DO iset=1,gto_basis_set%nset
00817        DO ishell=1,gto_basis_set%nshell(iset)
00818 
00819           l = gto_basis_set%l(ishell,iset)
00820 
00821           expzet = 0.5_dp*REAL(2*l + 3,dp)
00822 
00823           fnorm = 0.0_dp
00824 
00825           DO ipgf=1,gto_basis_set%npgf(iset)
00826              gcca = gto_basis_set%gcc(ipgf,ishell,iset)
00827              zeta = gto_basis_set%zet(ipgf,iset)
00828              DO jpgf=1,gto_basis_set%npgf(iset)
00829                 gccb = gto_basis_set%gcc(jpgf,ishell,iset)
00830                 zetb = gto_basis_set%zet(jpgf,iset)
00831                 fnorm = fnorm + gcca*gccb/(zeta + zetb)**expzet
00832              END DO
00833           END DO
00834 
00835           fnorm = 0.5_dp**l*pi**1.5_dp*fnorm
00836 
00837           DO icgf=gto_basis_set%first_cgf(ishell,iset),&
00838                gto_basis_set%last_cgf(ishell,iset)
00839              lx = gto_basis_set%lx(icgf)
00840              ly = gto_basis_set%ly(icgf)
00841              lz = gto_basis_set%lz(icgf)
00842              prefac = dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)
00843              gto_basis_set%norm_cgf(icgf) = 1.0_dp/SQRT(prefac*fnorm)
00844           END DO
00845 
00846        END DO
00847     END DO
00848 
00849   END SUBROUTINE init_norm_cgf_orb
00850 
00851 ! *****************************************************************************
00852   SUBROUTINE init_norm_cgf_orb_den(gto_basis_set,error)
00853 
00854     ! Initialise the normalization factors of the contracted Cartesian Gaussian
00855     ! functions used for frozen density representation.
00856 
00857     ! - Creation (21.09.2002,GT)
00858 
00859     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00860     TYPE(cp_error_type), INTENT(inout)       :: error
00861 
00862     INTEGER                                  :: icgf, ipgf, iset, ishell, l
00863     REAL(KIND=dp)                            :: expzet, gcca, prefac, zeta
00864 
00865 ! -------------------------------------------------------------------------
00866 
00867     DO iset=1,gto_basis_set%nset
00868        DO ishell=1,gto_basis_set%nshell(iset)
00869           l = gto_basis_set%l(ishell,iset)
00870           expzet = 0.5_dp*REAL(2*l + 3,dp)
00871           prefac = (1.0_dp/pi)**1.5_dp
00872           DO ipgf=1,gto_basis_set%npgf(iset)
00873              gcca = gto_basis_set%gcc(ipgf,ishell,iset)
00874              zeta = gto_basis_set%zet(ipgf,iset)
00875              gto_basis_set%gcc(ipgf,ishell,iset) = prefac*zeta**expzet*gcca
00876           END DO
00877           DO icgf=gto_basis_set%first_cgf(ishell,iset),&
00878                gto_basis_set%last_cgf(ishell,iset)
00879              gto_basis_set%norm_cgf(icgf) = 1.0_dp
00880           END DO
00881        END DO
00882     END DO
00883 
00884   END SUBROUTINE init_norm_cgf_orb_den
00885 
00886 ! *****************************************************************************
00887   SUBROUTINE init_orb_basis_set(gto_basis_set,error)
00888 
00889     ! Initialise a Gaussian-type orbital (GTO) basis set data set.
00890 
00891     ! - Creation (26.10.2000,MK)
00892 
00893     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00894     TYPE(cp_error_type), INTENT(inout)       :: error
00895 
00896     CHARACTER(len=*), PARAMETER :: routineN = 'init_orb_basis_set', 
00897       routineP = moduleN//':'//routineN
00898 
00899     INTEGER                                  :: handle
00900 
00901 ! -------------------------------------------------------------------------
00902 
00903     IF (.NOT.ASSOCIATED(gto_basis_set)) RETURN
00904 
00905     CALL timeset(routineN,handle)
00906 
00907     SELECT CASE (gto_basis_set%norm_type)
00908     CASE ( 0 )
00909        ! No normalisation requested
00910     CASE ( 1 )
00911        CALL init_norm_cgf_orb_den(gto_basis_set,error)
00912     CASE ( 2 )
00913        ! Normalise the primitive Gaussian functions
00914        CALL normalise_gcc_orb(gto_basis_set,error)
00915        ! Compute the normalization factors of the contracted Gaussian-type
00916        ! functions
00917        CALL init_norm_cgf_orb(gto_basis_set,error)
00918     CASE DEFAULT
00919        CALL stop_program(routineN,moduleN,__LINE__,&
00920             "Normalization method not specified")
00921     END SELECT
00922 
00923     ! Initialise the transformation matrices "pgf" -> "cgf"
00924 
00925     CALL init_cphi_and_sphi(gto_basis_set,error)
00926 
00927     CALL timestop(handle)
00928 
00929   END SUBROUTINE init_orb_basis_set
00930 
00931 ! *****************************************************************************
00932   SUBROUTINE normalise_gcc_orb(gto_basis_set,error)
00933 
00934     ! Normalise the primitive Cartesian Gaussian functions. The normalization
00935     ! factor is included in the Gaussian contraction coefficients.
00936 
00937     ! - Creation (20.08.1999,MK)
00938 
00939     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00940     TYPE(cp_error_type), INTENT(inout)       :: error
00941 
00942     INTEGER                                  :: ipgf, iset, ishell, l
00943     REAL(KIND=dp)                            :: expzet, gcca, prefac, zeta
00944 
00945 ! -------------------------------------------------------------------------
00946 
00947     DO iset=1,gto_basis_set%nset
00948        DO ishell=1,gto_basis_set%nshell(iset)
00949           l = gto_basis_set%l(ishell,iset)
00950           expzet = 0.25_dp*REAL(2*l + 3,dp)
00951           prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
00952           DO ipgf=1,gto_basis_set%npgf(iset)
00953              gcca = gto_basis_set%gcc(ipgf,ishell,iset)
00954              zeta = gto_basis_set%zet(ipgf,iset)
00955              gto_basis_set%gcc(ipgf,ishell,iset) = prefac*zeta**expzet*gcca
00956           END DO
00957        END DO
00958     END DO
00959 
00960   END SUBROUTINE normalise_gcc_orb
00961 
00962 ! *****************************************************************************
00963   SUBROUTINE read_gto_basis_set(element_symbol,basis_set_name,gto_basis_set,&
00964        para_env,dft_section,basis_section,error)
00965 
00966     ! Read a Gaussian-type orbital (GTO) basis set from the database file.
00967 
00968     ! - Creation (13.04.2000,MK)
00969 
00970     CHARACTER(LEN=*), INTENT(IN)             :: element_symbol, basis_set_name
00971     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
00972     TYPE(cp_para_env_type), POINTER          :: para_env
00973     TYPE(section_vals_type), POINTER         :: dft_section
00974     TYPE(section_vals_type), OPTIONAL, 
00975       POINTER                                :: basis_section
00976     TYPE(cp_error_type), INTENT(inout)       :: error
00977 
00978     CHARACTER(len=*), PARAMETER :: routineN = 'read_gto_basis_set', 
00979       routineP = moduleN//':'//routineN
00980 
00981     CHARACTER(len=20*default_string_length)  :: line_att
00982     CHARACTER(LEN=240)                       :: line
00983     CHARACTER(LEN=242)                       :: line2
00984     CHARACTER(len=default_path_length)       :: basis_set_file_name, tmp
00985     CHARACTER(LEN=default_string_length), 
00986       DIMENSION(:), POINTER                  :: cbasis
00987     CHARACTER(LEN=LEN(basis_set_name))       :: bsname
00988     CHARACTER(LEN=LEN(basis_set_name)+2)     :: bsname2
00989     CHARACTER(LEN=LEN(element_symbol))       :: symbol
00990     CHARACTER(LEN=LEN(element_symbol)+2)     :: symbol2
00991     INTEGER :: i, ibasis, ico, ipgf, irep, iset, ishell, istr, lshell, m, 
00992       maxco, maxl, maxpgf, maxshell, nbasis, ncgf, nmin, nset, nsgf, stat, 
00993       strlen1, strlen2
00994     INTEGER, DIMENSION(:), POINTER           :: lmax, lmin, npgf, nshell
00995     INTEGER, DIMENSION(:, :), POINTER        :: l, n
00996     LOGICAL                                  :: basis_found, failure, found, 
00997                                                 is_ok, match, read_from_input
00998     REAL(KIND=dp), DIMENSION(:, :), POINTER  :: zet
00999     REAL(KIND=dp), DIMENSION(:, :, :), 
01000       POINTER                                :: gcc
01001     TYPE(cp_parser_type), POINTER            :: parser
01002     TYPE(cp_sll_val_type), POINTER           :: list
01003     TYPE(val_type), POINTER                  :: val
01004 
01005     failure = .FALSE.
01006     line = ""
01007     line2 = ""
01008     symbol = ""
01009     symbol2 = ""
01010     bsname = ""
01011     bsname2 = ""
01012 
01013     nbasis = 1
01014 
01015     gto_basis_set%name = basis_set_name
01016     read_from_input = .FALSE.
01017     IF (PRESENT(basis_section)) THEN
01018        CALL section_vals_get(basis_section,explicit=read_from_input, error=error)
01019     END IF
01020     IF (.NOT.read_from_input) THEN
01021        CALL section_vals_val_get(dft_section,"BASIS_SET_FILE_NAME",&
01022             c_vals=cbasis,error=error)
01023        nbasis = SIZE(cbasis)
01024        DO ibasis = 1,nbasis
01025          NULLIFY(parser)
01026          basis_set_file_name = cbasis(ibasis)
01027          tmp=basis_set_file_name
01028          CALL uppercase(tmp)
01029          IF (INDEX(tmp,"MOLOPT").NE.0) CALL cite_reference(VandeVondele2007)
01030        END DO
01031     END IF
01032 
01033     ! Search for the requested basis set in the basis set file
01034     ! until the basis set is found or the end of file is reached
01035 
01036     basis_found = .FALSE.
01037     basis_loop:DO ibasis = 1,nbasis
01038       IF( basis_found ) EXIT basis_loop
01039       IF (.NOT.read_from_input) THEN
01040         NULLIFY(parser)
01041         basis_set_file_name = cbasis(ibasis)
01042         CALL parser_create(parser,basis_set_file_name,para_env=para_env,error=error)
01043       END IF
01044 
01045       bsname = basis_set_name
01046       symbol = element_symbol
01047       irep   = 0
01048 
01049       tmp=basis_set_name
01050       CALL uppercase(tmp)
01051       IF (INDEX(tmp,"MOLOPT").NE.0) CALL cite_reference(VandeVondele2007)
01052 
01053       nset = 0
01054       maxshell=0
01055       maxpgf=0
01056       maxco=0
01057       ncgf=0
01058       nsgf=0
01059       gto_basis_set%nset = nset
01060       gto_basis_set%ncgf = ncgf
01061       gto_basis_set%nsgf = nsgf
01062       CALL reallocate(gto_basis_set%lmax,1,nset)
01063       CALL reallocate(gto_basis_set%lmin,1,nset)
01064       CALL reallocate(gto_basis_set%npgf,1,nset)
01065       CALL reallocate(gto_basis_set%nshell,1,nset)
01066       CALL reallocate(gto_basis_set%n,1,maxshell,1,nset)
01067       CALL reallocate(gto_basis_set%l,1,maxshell,1,nset)
01068       CALL reallocate(gto_basis_set%zet,1,maxpgf,1,nset)
01069       CALL reallocate(gto_basis_set%gcc,1,maxpgf,1,maxshell,1,nset)
01070       CALL reallocate(gto_basis_set%set_radius,1,nset)
01071       CALL reallocate(gto_basis_set%pgf_radius,1,maxpgf,1,nset)
01072       CALL reallocate(gto_basis_set%first_cgf,1,maxshell,1,nset)
01073       CALL reallocate(gto_basis_set%first_sgf,1,maxshell,1,nset)
01074       CALL reallocate(gto_basis_set%last_cgf,1,maxshell,1,nset)
01075       CALL reallocate(gto_basis_set%last_sgf,1,maxshell,1,nset)
01076       CALL reallocate(gto_basis_set%ncgf_set,1,nset)
01077       CALL reallocate(gto_basis_set%nsgf_set,1,nset)
01078       CALL reallocate(gto_basis_set%cphi,1,maxco,1,ncgf)
01079       CALL reallocate(gto_basis_set%sphi,1,maxco,1,nsgf)
01080       CALL reallocate(gto_basis_set%lx,1,ncgf)
01081       CALL reallocate(gto_basis_set%ly,1,ncgf)
01082       CALL reallocate(gto_basis_set%lz,1,ncgf)
01083       CALL reallocate(gto_basis_set%m,1,nsgf)
01084       CALL reallocate(gto_basis_set%norm_cgf,1,ncgf)
01085 
01086       IF (tmp.NE."NONE") THEN
01087         search_loop: DO
01088 
01089            IF (read_from_input) THEN
01090               NULLIFY(list,val)
01091               found = .TRUE.
01092               CALL section_vals_list_get(basis_section,"_DEFAULT_KEYWORD_",list=list,error=error)
01093            ELSE
01094               CALL parser_search_string(parser,TRIM(bsname),.TRUE.,found,line,error=error)
01095            END IF
01096            IF (found) THEN
01097               CALL uppercase(symbol)
01098               CALL uppercase(bsname)
01099 
01100               IF (read_from_input) THEN
01101                  match = .TRUE.
01102               ELSE
01103                  match = .FALSE.
01104                  CALL uppercase(line)
01105                  ! Check both the element symbol and the basis set name
01106                  line2 = " "//line//" "
01107                  symbol2 = " "//TRIM(symbol)//" "
01108                  bsname2 = " "//TRIM(bsname)//" "
01109                  strlen1 = LEN_TRIM(symbol2) + 1
01110                  strlen2 = LEN_TRIM(bsname2) + 1
01111 
01112                  IF ( (INDEX(line2,symbol2(:strlen1)) > 0).AND.&
01113                       (INDEX(line2,bsname2(:strlen2)) > 0) ) match = .TRUE.
01114               END IF
01115               IF (match) THEN
01116                  NULLIFY (gcc,l,lmax,lmin,n,npgf,nshell,zet)
01117                  ! Read the basis set information
01118                  IF (read_from_input) THEN
01119                     is_ok=cp_sll_val_next(list,val,error=error)
01120                     IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,&
01121                          __LINE__,&
01122                          "Error reading the Basis set from input file!!")
01123                     CALL val_get(val,c_val=line_att,error=error)
01124                     READ(line_att,*)nset
01125                  ELSE
01126                     CALL parser_get_object(parser,nset,newline=.TRUE.,error=error)
01127                     IF (PRESENT(basis_section)) THEN
01128                        irep = irep + 1
01129                        WRITE(line_att,'(1X,I0)')nset
01130                        CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
01131                             c_val=TRIM(line_att), error=error)
01132                     END IF
01133                  END IF
01134 
01135                  CALL reallocate(npgf,1,nset)
01136                  CALL reallocate(nshell,1,nset)
01137                  CALL reallocate(lmax,1,nset)
01138                  CALL reallocate(lmin,1,nset)
01139                  CALL reallocate(n,1,1,1,nset)
01140 
01141                  maxl = 0
01142                  maxpgf = 0
01143                  maxshell = 0
01144 
01145                  DO iset=1,nset
01146                     IF (read_from_input) THEN
01147                        is_ok=cp_sll_val_next(list,val,error=error)
01148                        IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,&
01149                             __LINE__,&
01150                             "Error reading the Basis set from input file!!")
01151                        CALL val_get(val,c_val=line_att,error=error)
01152                        READ(line_att,*)n(1,iset)
01153                        CALL remove_word(line_att)
01154                        READ(line_att,*)lmin(iset)
01155                        CALL remove_word(line_att)
01156                        READ(line_att,*)lmax(iset)
01157                        CALL remove_word(line_att)
01158                        READ(line_att,*)npgf(iset)
01159                        CALL remove_word(line_att)
01160                     ELSE
01161                        line_att = ""
01162                        CALL parser_get_object(parser,n(1,iset),newline=.TRUE.,error=error)
01163                        CALL parser_get_object(parser,lmin(iset),error=error)
01164                        CALL parser_get_object(parser,lmax(iset),error=error)
01165                        CALL parser_get_object(parser,npgf(iset),error=error)
01166                        WRITE(line_att,'(4(1X,I0))')n(1,iset),lmin(iset),lmax(iset),npgf(iset)
01167                     END IF
01168                     maxl = MAX(maxl,lmax(iset))
01169                     IF (npgf(iset) > maxpgf) THEN
01170                        maxpgf = npgf(iset)
01171                        CALL reallocate(zet,1,maxpgf,1,nset)
01172                        CALL reallocate(gcc,1,maxpgf,1,maxshell,1,nset)
01173                     END IF
01174                     nshell(iset) = 0
01175                     DO lshell=lmin(iset),lmax(iset)
01176                        nmin = n(1,iset) + lshell - lmin(iset)
01177                        IF (read_from_input) THEN
01178                           READ(line_att,*)ishell
01179                           CALL remove_word(line_att)
01180                        ELSE
01181                           CALL parser_get_object(parser,ishell,error=error)
01182                           istr=LEN_TRIM(line_att)+1
01183                           WRITE(line_att(istr:),'(1X,I0)')ishell
01184                        END IF
01185                        nshell(iset) = nshell(iset) + ishell
01186                        IF (nshell(iset) > maxshell) THEN
01187                           maxshell = nshell(iset)
01188                           CALL reallocate(n,1,maxshell,1,nset)
01189                           CALL reallocate(l,1,maxshell,1,nset)
01190                           CALL reallocate(gcc,1,maxpgf,1,maxshell,1,nset)
01191                        END IF
01192                        DO i=1,ishell
01193                           n(nshell(iset)-ishell+i,iset) = nmin + i - 1
01194                           l(nshell(iset)-ishell+i,iset) = lshell
01195                        END DO
01196                     END DO
01197                     IF (.NOT.read_from_input) THEN
01198                        IF (PRESENT(basis_section)) THEN
01199                           irep = irep + 1
01200                           CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
01201                                c_val=TRIM(line_att), error=error)
01202                           line_att = ""
01203                        END IF
01204                     ELSE
01205                        IF (LEN_TRIM(line_att)/=0) CALL stop_program(routineN,moduleN,&
01206                             __LINE__,&
01207                             "Error reading the Basis from input file!!")
01208                     END IF
01209                     DO ipgf=1,npgf(iset)
01210                        IF (read_from_input) THEN
01211                           is_ok=cp_sll_val_next(list,val,error=error)
01212                           IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,&
01213                                __LINE__,&
01214                                "Error reading the Basis set from input file!!")
01215                           CALL val_get(val,c_val=line_att,error=error)
01216                           READ(line_att,*)zet(ipgf,iset),(gcc(ipgf,ishell,iset),ishell=1,nshell(iset))
01217                        ELSE
01218                           CALL parser_get_object(parser,zet(ipgf,iset),newline=.TRUE.,error=error)
01219                           DO ishell=1,nshell(iset)
01220                              CALL parser_get_object(parser,gcc(ipgf,ishell,iset),error=error)
01221                           END DO
01222                           IF (PRESENT(basis_section)) THEN
01223                              irep = irep + 1
01224                              WRITE(line_att,'(100E24.16)')zet(ipgf,iset),(gcc(ipgf,ishell,iset),ishell=1,nshell(iset))
01225                              CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
01226                                   c_val=TRIM(line_att), error=error)
01227                           END IF
01228                        END IF
01229                     END DO
01230                  END DO
01231 
01232                  ! Maximum angular momentum quantum number of the atomic kind
01233 
01234                  CALL init_orbital_pointers(maxl)
01235 
01236                  ! Allocate the global variables
01237 
01238                  gto_basis_set%nset = nset
01239                  CALL reallocate(gto_basis_set%lmax,1,nset)
01240                  CALL reallocate(gto_basis_set%lmin,1,nset)
01241                  CALL reallocate(gto_basis_set%npgf,1,nset)
01242                  CALL reallocate(gto_basis_set%nshell,1,nset)
01243                  CALL reallocate(gto_basis_set%n,1,maxshell,1,nset)
01244                  CALL reallocate(gto_basis_set%l,1,maxshell,1,nset)
01245                  CALL reallocate(gto_basis_set%zet,1,maxpgf,1,nset)
01246                  CALL reallocate(gto_basis_set%gcc,1,maxpgf,1,maxshell,1,nset)
01247 
01248                  ! Copy the basis set information into the data structure
01249 
01250                  DO iset=1,nset
01251                     gto_basis_set%lmax(iset) = lmax(iset)
01252                     gto_basis_set%lmin(iset) = lmin(iset)
01253                     gto_basis_set%npgf(iset) = npgf(iset)
01254                     gto_basis_set%nshell(iset) = nshell(iset)
01255                     DO ishell=1,nshell(iset)
01256                        gto_basis_set%n(ishell,iset) = n(ishell,iset)
01257                        gto_basis_set%l(ishell,iset) = l(ishell,iset)
01258                        DO ipgf=1,npgf(iset)
01259                           gto_basis_set%gcc(ipgf,ishell,iset) = gcc(ipgf,ishell,iset)
01260                        END DO
01261                     END DO
01262                     DO ipgf=1,npgf(iset)
01263                        gto_basis_set%zet(ipgf,iset) = zet(ipgf,iset)
01264                     END DO
01265                  END DO
01266 
01267                  ! Initialise the depending atomic kind information
01268 
01269                  CALL reallocate(gto_basis_set%set_radius,1,nset)
01270                  CALL reallocate(gto_basis_set%pgf_radius,1,maxpgf,1,nset)
01271                  CALL reallocate(gto_basis_set%first_cgf,1,maxshell,1,nset)
01272                  CALL reallocate(gto_basis_set%first_sgf,1,maxshell,1,nset)
01273                  CALL reallocate(gto_basis_set%last_cgf,1,maxshell,1,nset)
01274                  CALL reallocate(gto_basis_set%last_sgf,1,maxshell,1,nset)
01275                  CALL reallocate(gto_basis_set%ncgf_set,1,nset)
01276                  CALL reallocate(gto_basis_set%nsgf_set,1,nset)
01277 
01278                  maxco = 0
01279                  ncgf = 0
01280                  nsgf = 0
01281 
01282                  DO iset=1,nset
01283                     gto_basis_set%ncgf_set(iset) = 0
01284                     gto_basis_set%nsgf_set(iset) = 0
01285                     DO ishell=1,nshell(iset)
01286                        lshell = gto_basis_set%l(ishell,iset)
01287                        gto_basis_set%first_cgf(ishell,iset) = ncgf + 1
01288                        ncgf = ncgf + nco(lshell)
01289                        gto_basis_set%last_cgf(ishell,iset) = ncgf
01290                        gto_basis_set%ncgf_set(iset) =&
01291                             gto_basis_set%ncgf_set(iset) + nco(lshell)
01292                        gto_basis_set%first_sgf(ishell,iset) = nsgf + 1
01293                        nsgf = nsgf + nso(lshell)
01294                        gto_basis_set%last_sgf(ishell,iset) = nsgf
01295                        gto_basis_set%nsgf_set(iset) =&
01296                             gto_basis_set%nsgf_set(iset) + nso(lshell)
01297                     END DO
01298                     maxco = MAX(maxco,npgf(iset)*ncoset(lmax(iset)))
01299                  END DO
01300 
01301                  gto_basis_set%ncgf = ncgf
01302                  gto_basis_set%nsgf = nsgf
01303 
01304                  CALL reallocate(gto_basis_set%cphi,1,maxco,1,ncgf)
01305                  CALL reallocate(gto_basis_set%sphi,1,maxco,1,nsgf)
01306                  CALL reallocate(gto_basis_set%lx,1,ncgf)
01307                  CALL reallocate(gto_basis_set%ly,1,ncgf)
01308                  CALL reallocate(gto_basis_set%lz,1,ncgf)
01309                  CALL reallocate(gto_basis_set%m,1,nsgf)
01310                  CALL reallocate(gto_basis_set%norm_cgf,1,ncgf)
01311                  ALLOCATE (gto_basis_set%cgf_symbol(ncgf),STAT=stat)
01312                  CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01313 
01314                  ALLOCATE (gto_basis_set%sgf_symbol(nsgf),STAT=stat)
01315                  CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01316 
01317                  ncgf = 0
01318                  nsgf = 0
01319 
01320                  DO iset=1,nset
01321                     DO ishell=1,nshell(iset)
01322                        lshell = gto_basis_set%l(ishell,iset)
01323                        DO ico=ncoset(lshell-1)+1,ncoset(lshell)
01324                           ncgf = ncgf + 1
01325                           gto_basis_set%lx(ncgf) = indco(1,ico)
01326                           gto_basis_set%ly(ncgf) = indco(2,ico)
01327                           gto_basis_set%lz(ncgf) = indco(3,ico)
01328                           gto_basis_set%cgf_symbol(ncgf) =&
01329                                cgf_symbol(n(ishell,iset),(/gto_basis_set%lx(ncgf),&
01330                                gto_basis_set%ly(ncgf),&
01331                                gto_basis_set%lz(ncgf)/))
01332                        END DO
01333                        DO m=-lshell,lshell
01334                           nsgf = nsgf + 1
01335                           gto_basis_set%m(nsgf) = m
01336                           gto_basis_set%sgf_symbol(nsgf) =&
01337                                sgf_symbol(n(ishell,iset),lshell,m)
01338                        END DO
01339                     END DO
01340                  END DO
01341 
01342                  DEALLOCATE (gcc,l,lmax,lmin,n,npgf,nshell,zet,STAT=stat)
01343                  CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01344 
01345                  basis_found = .TRUE.
01346                  EXIT search_loop
01347               END IF
01348            ELSE
01349              EXIT search_loop
01350            END IF
01351         END DO search_loop
01352       ELSE
01353         match=.FALSE.
01354       ENDIF
01355 
01356       IF (.NOT.read_from_input .AND. basis_found) THEN
01357         CALL parser_release(parser,error=error)
01358         IF ((match).AND.(PRESENT(basis_section))) THEN
01359            ! Dump the read basis set in the basis section
01360            irep = irep + 1
01361            WRITE(line_att,'(A)')"         # Basis set name: "//bsname2(:strlen2)//" for symbol: "//symbol2(:strlen1)
01362            CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
01363                 c_val=TRIM(line_att), error=error)
01364            irep = irep + 1
01365           WRITE(line_att,'(A)')"         # Basis set read from the basis set filename: "//TRIM(basis_set_file_name)
01366             CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
01367                c_val=TRIM(line_att), error=error)
01368         END IF
01369       END IF
01370       IF ( .NOT. basis_found ) THEN
01371         CALL parser_release(parser,error=error)
01372       END IF
01373 
01374     END DO basis_loop
01375 
01376     IF (.NOT. read_from_input .AND. ( tmp .NE. "NONE") ) THEN
01377       IF( .NOT. basis_found ) THEN
01378         basis_set_file_name = ""
01379         DO ibasis = 1,nbasis
01380           basis_set_file_name = TRIM(basis_set_file_name)//"<"//TRIM(cbasis(ibasis))//"> "
01381         END DO
01382         CALL stop_program(routineN,moduleN,__LINE__,&
01383              "The requested basis set <"//TRIM(bsname)//&
01384              "> for element <"//TRIM(symbol)//"> was not "//&
01385              "found in the basis set files "//&
01386              TRIM(basis_set_file_name),para_env)
01387       END IF
01388     END IF
01389 
01390   END SUBROUTINE read_gto_basis_set
01391 
01392 ! *****************************************************************************
01393   SUBROUTINE set_gto_basis_set(gto_basis_set,name,norm_type,kind_radius,ncgf,&
01394        nset,nsgf,cgf_symbol,sgf_symbol,norm_cgf,set_radius,lmax,&
01395        lmin,lx,ly,lz,m,ncgf_set,npgf,nsgf_set,nshell,&
01396        cphi,pgf_radius,sphi,zet,first_cgf,first_sgf,l,&
01397        last_cgf,last_sgf,n,gcc,short_kind_radius)
01398 
01399     ! Set the components of Gaussian-type orbital (GTO) basis set data set.
01400 
01401     ! - Creation (10.01.2002,MK)
01402 
01403     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
01404     CHARACTER(LEN=default_string_length), 
01405       INTENT(IN), OPTIONAL                   :: name
01406     INTEGER, INTENT(IN), OPTIONAL            :: norm_type
01407     REAL(KIND=dp), INTENT(IN), OPTIONAL      :: kind_radius
01408     INTEGER, INTENT(IN), OPTIONAL            :: ncgf, nset, nsgf
01409     CHARACTER(LEN=12), DIMENSION(:), 
01410       OPTIONAL, POINTER                      :: cgf_symbol
01411     CHARACTER(LEN=6), DIMENSION(:), 
01412       OPTIONAL, POINTER                      :: sgf_symbol
01413     REAL(KIND=dp), DIMENSION(:), OPTIONAL, 
01414       POINTER                                :: norm_cgf, set_radius
01415     INTEGER, DIMENSION(:), OPTIONAL, POINTER :: lmax, lmin, lx, ly, lz, m, 
01416                                                 ncgf_set, npgf, nsgf_set, 
01417                                                 nshell
01418     REAL(KIND=dp), DIMENSION(:, :), 
01419       OPTIONAL, POINTER                      :: cphi, pgf_radius, sphi, zet
01420     INTEGER, DIMENSION(:, :), OPTIONAL, 
01421       POINTER                                :: first_cgf, first_sgf, l, 
01422                                                 last_cgf, last_sgf, n
01423     REAL(KIND=dp), DIMENSION(:, :, :), 
01424       OPTIONAL, POINTER                      :: gcc
01425     REAL(KIND=dp), INTENT(IN), OPTIONAL      :: short_kind_radius
01426 
01427     CHARACTER(len=*), PARAMETER :: routineN = 'set_gto_basis_set', 
01428       routineP = moduleN//':'//routineN
01429 
01430 ! -------------------------------------------------------------------------
01431 
01432     IF (ASSOCIATED(gto_basis_set)) THEN
01433        IF (PRESENT(name)) gto_basis_set%name = name
01434        IF (PRESENT(norm_type)) gto_basis_set%norm_type = norm_type
01435        IF (PRESENT(kind_radius)) gto_basis_set%kind_radius = kind_radius
01436        IF (PRESENT(short_kind_radius)) gto_basis_set%short_kind_radius = short_kind_radius
01437        IF (PRESENT(ncgf)) gto_basis_set%ncgf = ncgf
01438        IF (PRESENT(nset)) gto_basis_set%nset = nset
01439        IF (PRESENT(nsgf)) gto_basis_set%nsgf = nsgf
01440        IF (PRESENT(cgf_symbol)) gto_basis_set%cgf_symbol(:) = cgf_symbol(:)
01441        IF (PRESENT(sgf_symbol)) gto_basis_set%sgf_symbol(:) = sgf_symbol(:)
01442        IF (PRESENT(norm_cgf)) gto_basis_set%norm_cgf(:) = norm_cgf(:)
01443        IF (PRESENT(set_radius)) gto_basis_set%set_radius(:) = set_radius(:)
01444        IF (PRESENT(lmax)) gto_basis_set%lmax(:) = lmax(:)
01445        IF (PRESENT(lmin)) gto_basis_set%lmin(:) = lmin(:)
01446        IF (PRESENT(lx)) gto_basis_set%lx(:) = lx(:)
01447        IF (PRESENT(ly)) gto_basis_set%ly(:) = ly(:)
01448        IF (PRESENT(lz)) gto_basis_set%lz(:) = lz(:)
01449        IF (PRESENT(m)) gto_basis_set%m(:) = m(:)
01450        IF (PRESENT(ncgf_set)) gto_basis_set%ncgf_set(:) = ncgf_set(:)
01451        IF (PRESENT(npgf)) gto_basis_set%npgf(:) = npgf(:)
01452        IF (PRESENT(nsgf_set)) gto_basis_set%nsgf_set(:) = nsgf_set(:)
01453        IF (PRESENT(nshell)) gto_basis_set%nshell(:) = nshell(:)
01454        IF (PRESENT(cphi)) gto_basis_set%cphi(:,:) = cphi(:,:)
01455        IF (PRESENT(pgf_radius)) gto_basis_set%pgf_radius(:,:) = pgf_radius(:,:)
01456        IF (PRESENT(sphi)) gto_basis_set%sphi(:,:) = sphi(:,:)
01457        IF (PRESENT(zet)) gto_basis_set%zet(:,:) = zet(:,:)
01458        IF (PRESENT(first_cgf)) gto_basis_set%first_cgf(:,:) = first_cgf(:,:)
01459        IF (PRESENT(first_sgf)) gto_basis_set%first_sgf(:,:) = first_sgf(:,:)
01460        IF (PRESENT(l)) l(:,:) = gto_basis_set%l(:,:)
01461        IF (PRESENT(last_cgf)) gto_basis_set%last_cgf(:,:) = last_cgf(:,:)
01462        IF (PRESENT(last_sgf)) gto_basis_set%last_sgf(:,:) = last_sgf(:,:)
01463        IF (PRESENT(n)) gto_basis_set%n(:,:) = n(:,:)
01464        IF (PRESENT(gcc)) gto_basis_set%gcc(:,:,:) = gcc(:,:,:)
01465     ELSE
01466        CALL stop_program(routineN,moduleN,__LINE__,&
01467             "The pointer gto_basis_set is not associated")
01468     END IF
01469 
01470   END SUBROUTINE set_gto_basis_set
01471 
01472 ! *****************************************************************************
01473   SUBROUTINE write_gto_basis_set(gto_basis_set,output_unit,header,error)
01474 
01475     ! Write a Gaussian-type orbital (GTO) basis set data set to the output
01476     ! unit.
01477 
01478     ! - Creation (09.01.2002,MK)
01479 
01480     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
01481     INTEGER, INTENT(in)                      :: output_unit
01482     CHARACTER(len=*), OPTIONAL               :: header
01483     TYPE(cp_error_type), INTENT(inout)       :: error
01484 
01485     INTEGER                                  :: ipgf, iset, ishell
01486 
01487 ! -------------------------------------------------------------------------
01488 
01489     IF (ASSOCIATED(gto_basis_set).AND.(output_unit>0)) THEN
01490 
01491        IF(PRESENT(header)) THEN
01492          WRITE (UNIT=output_unit,FMT="(/,T6,A,T41,A40)")&
01493               TRIM(header),TRIM(gto_basis_set%name)
01494        END IF
01495 
01496        WRITE (UNIT=output_unit,FMT="(/,(T8,A,T71,I10))")&
01497             "Number of orbital shell sets:            ",&
01498             gto_basis_set%nset,&
01499             "Number of orbital shells:                ",&
01500             SUM(gto_basis_set%nshell(:)),&
01501             "Number of primitive Cartesian functions: ",&
01502             SUM(gto_basis_set%npgf(:)),&
01503             "Number of Cartesian basis functions:     ",&
01504             gto_basis_set%ncgf,&
01505             "Number of spherical basis functions:     ",&
01506             gto_basis_set%nsgf,&
01507             "Norm type:                               ",&
01508             gto_basis_set%norm_type
01509 
01510        WRITE (UNIT=output_unit,FMT="(/,T6,A,T41,A40,/,/,T25,A)")&
01511             "GTO basis set information for",TRIM(gto_basis_set%name),&
01512             "Set   Shell     n   l            Exponent    Coefficient"
01513 
01514        DO iset=1,gto_basis_set%nset
01515           WRITE (UNIT=output_unit,FMT="(A)") ""
01516           DO ishell=1,gto_basis_set%nshell(iset)
01517              WRITE (UNIT=output_unit,&
01518                   FMT="(T25,I3,4X,I4,4X,I2,2X,I2,(T51,2F15.6))")&
01519                   iset,ishell,&
01520                   gto_basis_set%n(ishell,iset),&
01521                   gto_basis_set%l(ishell,iset),&
01522                   (gto_basis_set%zet(ipgf,iset),&
01523                   gto_basis_set%gcc(ipgf,ishell,iset),&
01524                   ipgf=1,gto_basis_set%npgf(iset))
01525           END DO
01526        END DO
01527 
01528     END IF
01529 
01530   END SUBROUTINE write_gto_basis_set
01531 
01532 ! *****************************************************************************
01533 
01534   SUBROUTINE write_orb_basis_set(orb_basis_set,output_unit,header,error)
01535 
01536     ! Write a Gaussian-type orbital (GTO) basis set data set to the output
01537     ! unit.
01538 
01539     ! - Creation (09.01.2002,MK)
01540 
01541     TYPE(gto_basis_set_type), POINTER        :: orb_basis_set
01542     INTEGER, INTENT(in)                      :: output_unit
01543     CHARACTER(len=*), OPTIONAL               :: header
01544     TYPE(cp_error_type), INTENT(inout)       :: error
01545 
01546     INTEGER                                  :: icgf, ico, ipgf, iset, ishell
01547 
01548     IF (ASSOCIATED(orb_basis_set).AND.(output_unit>0)) THEN
01549        IF(PRESENT(header)) THEN
01550          WRITE (UNIT=output_unit,FMT="(/,T6,A,T41,A40)")&
01551               TRIM(header),TRIM(orb_basis_set%name)
01552        END IF
01553 
01554        WRITE (UNIT=output_unit,FMT="(/,(T8,A,T71,I10))")&
01555             "Number of orbital shell sets:            ",&
01556             orb_basis_set%nset,&
01557             "Number of orbital shells:                ",&
01558             SUM(orb_basis_set%nshell(:)),&
01559             "Number of primitive Cartesian functions: ",&
01560             SUM(orb_basis_set%npgf(:)),&
01561             "Number of Cartesian basis functions:     ",&
01562             orb_basis_set%ncgf,&
01563             "Number of spherical basis functions:     ",&
01564             orb_basis_set%nsgf,&
01565             "Norm type:                               ",&
01566             orb_basis_set%norm_type
01567 
01568        WRITE (UNIT=output_unit,FMT="(/,T8,A,/,/,T25,A)")&
01569             "Normalised Cartesian orbitals:",&
01570             "Set   Shell   Orbital            Exponent    Coefficient"
01571 
01572        icgf = 0
01573 
01574        DO iset=1,orb_basis_set%nset
01575           DO ishell=1,orb_basis_set%nshell(iset)
01576              WRITE (UNIT=output_unit,FMT="(A)") ""
01577              DO ico=1,nco(orb_basis_set%l(ishell,iset))
01578                 icgf = icgf + 1
01579                 WRITE (UNIT=output_unit,&
01580                      FMT="(T25,I3,4X,I4,3X,A12,(T51,2F15.6))")&
01581                      iset,ishell,orb_basis_set%cgf_symbol(icgf),&
01582                      (orb_basis_set%zet(ipgf,iset),&
01583                      orb_basis_set%norm_cgf(icgf)*&
01584                      orb_basis_set%gcc(ipgf,ishell,iset),&
01585                      ipgf=1,orb_basis_set%npgf(iset))
01586              END DO
01587           END DO
01588        END DO
01589     END IF
01590 
01591   END SUBROUTINE write_orb_basis_set
01592 
01593 ! *****************************************************************************
01594   SUBROUTINE allocate_sto_basis_set(sto_basis_set, error)
01595 
01596     TYPE(sto_basis_set_type), POINTER        :: sto_basis_set
01597     TYPE(cp_error_type), INTENT(inout)       :: error
01598 
01599     CHARACTER(len=*), PARAMETER :: routineN = 'allocate_sto_basis_set', 
01600       routineP = moduleN//':'//routineN
01601 
01602     INTEGER                                  :: stat
01603     LOGICAL                                  :: failure
01604 
01605 ! -------------------------------------------------------------------------
01606 
01607     failure = .FALSE.
01608     CALL deallocate_sto_basis_set(sto_basis_set, error)
01609 
01610     ALLOCATE (sto_basis_set,STAT=stat)
01611     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01612 
01613     sto_basis_set%name="NONAME"
01614     NULLIFY (sto_basis_set%symbol)
01615     NULLIFY (sto_basis_set%nq)
01616     NULLIFY (sto_basis_set%lq)
01617     NULLIFY (sto_basis_set%zet)
01618 
01619   END SUBROUTINE allocate_sto_basis_set
01620 
01621 ! *****************************************************************************
01622   SUBROUTINE deallocate_sto_basis_set(sto_basis_set, error)
01623 
01624     TYPE(sto_basis_set_type), POINTER        :: sto_basis_set
01625     TYPE(cp_error_type), INTENT(inout)       :: error
01626 
01627     CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_sto_basis_set', 
01628       routineP = moduleN//':'//routineN
01629 
01630     INTEGER                                  :: stat
01631     LOGICAL                                  :: failure
01632 
01633 ! -------------------------------------------------------------------------
01634 
01635     failure = .FALSE.
01636     IF (ASSOCIATED(sto_basis_set)) THEN
01637        IF (ASSOCIATED(sto_basis_set%symbol)) THEN
01638           DEALLOCATE (sto_basis_set%symbol,STAT=stat)
01639           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01640        END IF
01641        IF (ASSOCIATED(sto_basis_set%nq)) THEN
01642           DEALLOCATE (sto_basis_set%nq,STAT=stat)
01643           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01644        END IF
01645        IF (ASSOCIATED(sto_basis_set%lq)) THEN
01646           DEALLOCATE (sto_basis_set%lq,STAT=stat)
01647           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01648        END IF
01649        IF (ASSOCIATED(sto_basis_set%zet)) THEN
01650           DEALLOCATE (sto_basis_set%zet,STAT=stat)
01651           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01652        END IF
01653 
01654        DEALLOCATE (sto_basis_set,STAT=stat)
01655        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01656     END IF
01657   END SUBROUTINE deallocate_sto_basis_set
01658 
01659 ! *****************************************************************************
01660   SUBROUTINE get_sto_basis_set(sto_basis_set,name,nshell,symbol,nq,lq,zet,maxlq,numsto)
01661 
01662     TYPE(sto_basis_set_type), POINTER        :: sto_basis_set
01663     CHARACTER(LEN=default_string_length), 
01664       INTENT(OUT), OPTIONAL                  :: name
01665     INTEGER, INTENT(OUT), OPTIONAL           :: nshell
01666     CHARACTER(LEN=6), DIMENSION(:), 
01667       OPTIONAL, POINTER                      :: symbol
01668     INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nq, lq
01669     REAL(KIND=dp), DIMENSION(:), OPTIONAL, 
01670       POINTER                                :: zet
01671     INTEGER, INTENT(OUT), OPTIONAL           :: maxlq, numsto
01672 
01673     CHARACTER(len=*), PARAMETER :: routineN = 'get_sto_basis_set', 
01674       routineP = moduleN//':'//routineN
01675 
01676     INTEGER                                  :: iset
01677 
01678 ! -------------------------------------------------------------------------
01679 
01680     IF (ASSOCIATED(sto_basis_set)) THEN
01681        IF (PRESENT(name)) name = sto_basis_set%name
01682        IF (PRESENT(nshell)) nshell = sto_basis_set%nshell
01683        IF (PRESENT(symbol)) symbol => sto_basis_set%symbol
01684        IF (PRESENT(nq)) nq => sto_basis_set%nq
01685        IF (PRESENT(lq)) lq => sto_basis_set%lq
01686        IF (PRESENT(zet)) zet => sto_basis_set%zet
01687        IF (PRESENT(maxlq)) THEN
01688           maxlq = MAXVAL ( sto_basis_set%lq (1:sto_basis_set%nshell) )
01689        END IF
01690        IF (PRESENT(numsto)) THEN
01691           numsto = 0
01692           DO iset=1,sto_basis_set%nshell
01693              numsto = numsto + 2*sto_basis_set%lq(iset)+1
01694           END DO
01695        END IF
01696     ELSE
01697        CALL stop_program(routineN,moduleN,__LINE__,&
01698             "The pointer sto_basis_set is not associated")
01699     END IF
01700 
01701   END SUBROUTINE get_sto_basis_set
01702 
01703 ! *****************************************************************************
01704   SUBROUTINE set_sto_basis_set(sto_basis_set,name,nshell,symbol,nq,lq,zet)
01705 
01706     TYPE(sto_basis_set_type), POINTER        :: sto_basis_set
01707     CHARACTER(LEN=default_string_length), 
01708       INTENT(IN), OPTIONAL                   :: name
01709     INTEGER, INTENT(IN), OPTIONAL            :: nshell
01710     CHARACTER(LEN=6), DIMENSION(:), 
01711       OPTIONAL, POINTER                      :: symbol
01712     INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nq, lq
01713     REAL(KIND=dp), DIMENSION(:), OPTIONAL, 
01714       POINTER                                :: zet
01715 
01716     CHARACTER(len=*), PARAMETER :: routineN = 'set_sto_basis_set', 
01717       routineP = moduleN//':'//routineN
01718 
01719     INTEGER                                  :: ns
01720 
01721 ! -------------------------------------------------------------------------
01722 
01723     IF (ASSOCIATED(sto_basis_set)) THEN
01724        IF (PRESENT(name)) sto_basis_set%name = name
01725        IF (PRESENT(nshell)) sto_basis_set%nshell = nshell
01726        IF (PRESENT(symbol)) THEN
01727           ns = SIZE(symbol)
01728           IF (ASSOCIATED(sto_basis_set%symbol)) DEALLOCATE (sto_basis_set%symbol)
01729           ALLOCATE (sto_basis_set%symbol(1:ns))
01730           sto_basis_set%symbol(:) = symbol(:)
01731        END IF
01732        IF (PRESENT(nq)) THEN
01733           ns = SIZE(nq)
01734           CALL reallocate(sto_basis_set%nq,1,ns)
01735           sto_basis_set%nq = nq(:)
01736        END IF
01737        IF (PRESENT(lq)) THEN
01738           ns = SIZE(lq)
01739           CALL reallocate(sto_basis_set%lq,1,ns)
01740           sto_basis_set%lq = lq(:)
01741        END IF
01742        IF (PRESENT(zet)) THEN
01743           ns = SIZE(zet)
01744           CALL reallocate(sto_basis_set%zet,1,ns)
01745           sto_basis_set%zet = zet(:)
01746        END IF
01747     ELSE
01748        CALL stop_program(routineN,moduleN,__LINE__,&
01749             "The pointer sto_basis_set is not associated")
01750     END IF
01751 
01752   END SUBROUTINE set_sto_basis_set
01753 
01754 ! *****************************************************************************
01755   SUBROUTINE create_gto_from_sto_basis(sto_basis_set,gto_basis_set,ngauss,error)
01756 
01757     TYPE(sto_basis_set_type), POINTER        :: sto_basis_set
01758     TYPE(gto_basis_set_type), POINTER        :: gto_basis_set
01759     INTEGER, OPTIONAL                        :: ngauss
01760     TYPE(cp_error_type), INTENT(inout)       :: error
01761 
01762     CHARACTER(len=*), PARAMETER :: routineN = 'create_gto_from_sto_basis', 
01763       routineP = moduleN//':'//routineN
01764     INTEGER, PARAMETER                       :: maxng = 6
01765 
01766     CHARACTER(LEN=default_string_length)     :: name, sng
01767     INTEGER                                  :: ico, ipgf, iset, lshell, m, 
01768                                                 maxco, maxl, ncgf, ng = 6, 
01769                                                 np, nset, nsgf, nshell, stat
01770     INTEGER, DIMENSION(:), POINTER           :: lq, nq
01771     LOGICAL                                  :: failure
01772     REAL(KIND=dp), DIMENSION(:), POINTER     :: zet
01773     REAL(KIND=dp), DIMENSION(maxng)          :: gcc, zetg
01774 
01775     failure = .FALSE.
01776     IF (PRESENT(ngauss)) ng=ngauss
01777     IF (ng > maxng) CALL stop_program(routineN,moduleN,__LINE__,&
01778          "Too many Gaussian primitives requested")
01779 
01780     CALL allocate_gto_basis_set(gto_basis_set,error)
01781 
01782     CALL get_sto_basis_set(sto_basis_set,name=name,nshell=nshell,nq=nq,&
01783          lq=lq,zet=zet)
01784 
01785     maxl = MAXVAL(lq)
01786     CALL init_orbital_pointers(maxl)
01787 
01788     CALL integer_to_string(ng,sng)
01789     gto_basis_set%name = TRIM(name)//"_STO-"//TRIM(sng)//"G"
01790 
01791     nset = nshell
01792     gto_basis_set%nset = nset
01793     CALL reallocate(gto_basis_set%lmax,1,nset)
01794     CALL reallocate(gto_basis_set%lmin,1,nset)
01795     CALL reallocate(gto_basis_set%npgf,1,nset)
01796     CALL reallocate(gto_basis_set%nshell,1,nset)
01797     CALL reallocate(gto_basis_set%n,1,1,1,nset)
01798     CALL reallocate(gto_basis_set%l,1,1,1,nset)
01799     CALL reallocate(gto_basis_set%zet,1,ng,1,nset)
01800     CALL reallocate(gto_basis_set%gcc,1,ng,1,1,1,nset)
01801 
01802     DO iset=1,nset
01803        CALL get_sto_ng ( zet(iset), ng, nq(iset), lq(iset), zetg, gcc, error )
01804        gto_basis_set%lmax(iset) = lq(iset)
01805        gto_basis_set%lmin(iset) = lq(iset)
01806        gto_basis_set%npgf(iset) = ng
01807        gto_basis_set%nshell(iset) = 1
01808        gto_basis_set%n(1,iset) = lq(iset)+1
01809        gto_basis_set%l(1,iset) = lq(iset)
01810        DO ipgf=1,ng
01811           gto_basis_set%gcc(ipgf,1,iset) = gcc(ipgf)
01812           gto_basis_set%zet(ipgf,iset) = zetg(ipgf)
01813        END DO
01814     END DO
01815 
01816     CALL reallocate(gto_basis_set%set_radius,1,nset)
01817     CALL reallocate(gto_basis_set%pgf_radius,1,ng,1,nset)
01818     CALL reallocate(gto_basis_set%first_cgf,1,1,1,nset)
01819     CALL reallocate(gto_basis_set%first_sgf,1,1,1,nset)
01820     CALL reallocate(gto_basis_set%last_cgf,1,1,1,nset)
01821     CALL reallocate(gto_basis_set%last_sgf,1,1,1,nset)
01822     CALL reallocate(gto_basis_set%ncgf_set,1,nset)
01823     CALL reallocate(gto_basis_set%nsgf_set,1,nset)
01824 
01825     maxco = 0
01826     ncgf = 0
01827     nsgf = 0
01828 
01829     DO iset=1,nset
01830        gto_basis_set%ncgf_set(iset) = 0
01831        gto_basis_set%nsgf_set(iset) = 0
01832        lshell = gto_basis_set%l(1,iset)
01833        gto_basis_set%first_cgf(1,iset) = ncgf + 1
01834        ncgf = ncgf + nco(lshell)
01835        gto_basis_set%last_cgf(1,iset) = ncgf
01836        gto_basis_set%ncgf_set(iset) =&
01837             gto_basis_set%ncgf_set(iset) + nco(lshell)
01838        gto_basis_set%first_sgf(1,iset) = nsgf + 1
01839        nsgf = nsgf + nso(lshell)
01840        gto_basis_set%last_sgf(1,iset) = nsgf
01841        gto_basis_set%nsgf_set(iset) =&
01842             gto_basis_set%nsgf_set(iset) + nso(lshell)
01843        maxco = MAX(maxco,ng*ncoset(lshell))
01844     END DO
01845 
01846     gto_basis_set%ncgf = ncgf
01847     gto_basis_set%nsgf = nsgf
01848 
01849     CALL reallocate(gto_basis_set%cphi,1,maxco,1,ncgf)
01850     CALL reallocate(gto_basis_set%sphi,1,maxco,1,nsgf)
01851     CALL reallocate(gto_basis_set%lx,1,ncgf)
01852     CALL reallocate(gto_basis_set%ly,1,ncgf)
01853     CALL reallocate(gto_basis_set%lz,1,ncgf)
01854     CALL reallocate(gto_basis_set%m,1,nsgf)
01855     CALL reallocate(gto_basis_set%norm_cgf,1,ncgf)
01856     ALLOCATE (gto_basis_set%cgf_symbol(ncgf),STAT=stat)
01857     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01858     ALLOCATE (gto_basis_set%sgf_symbol(nsgf),STAT=stat)
01859     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01860 
01861     ncgf = 0
01862     nsgf = 0
01863 
01864     DO iset=1,nset
01865        lshell = gto_basis_set%l(1,iset)
01866        np = lshell + 1
01867        DO ico=ncoset(lshell-1)+1,ncoset(lshell)
01868           ncgf = ncgf + 1
01869           gto_basis_set%lx(ncgf) = indco(1,ico)
01870           gto_basis_set%ly(ncgf) = indco(2,ico)
01871           gto_basis_set%lz(ncgf) = indco(3,ico)
01872           gto_basis_set%cgf_symbol(ncgf) =&
01873                cgf_symbol(np,(/gto_basis_set%lx(ncgf),&
01874                gto_basis_set%ly(ncgf),&
01875                gto_basis_set%lz(ncgf)/))
01876        END DO
01877        DO m=-lshell,lshell
01878           nsgf = nsgf + 1
01879           gto_basis_set%m(nsgf) = m
01880           gto_basis_set%sgf_symbol(nsgf) = sgf_symbol(np,lshell,m)
01881        END DO
01882     END DO
01883 
01884     gto_basis_set%norm_type = -1
01885 
01886   END SUBROUTINE create_gto_from_sto_basis
01887 
01888  ! *****************************************************************************
01889   FUNCTION srules(z,ne,n,l)
01890     ! Slater rules
01891     INTEGER                                  :: z
01892     INTEGER, DIMENSION(4, 7)                 :: ne
01893     INTEGER                                  :: n, l
01894     REAL(dp)                                 :: srules
01895 
01896     REAL(dp), DIMENSION(7), PARAMETER :: 
01897       xns = (/ 1.0_dp,2.0_dp,3.0_dp,3.7_dp,4.0_dp,4.2_dp,4.4_dp /)
01898 
01899     INTEGER                                  :: i, l1, l2, m, m1, m2
01900     REAL(dp)                                 :: s
01901 
01902     s = 0.0_dp
01903     ! The complete shell
01904     l1=l+1
01905     IF(l1 == 1) l2=2
01906     IF(l1 == 2) l2=1
01907     IF(l1 == 3) l2=4
01908     IF(l1 == 4) l2=3
01909     ! Rule a) no contribution from shells further out
01910     ! Rule b) 0.35 (1s 0.3) from each other electron in the same shell
01911     IF(n == 1) THEN
01912        m=ne(1,1)
01913        s=s+0.3_dp*REAL(m-1,dp)
01914     ELSE
01915        m=ne(l1,n)+ne(l2,n)
01916        s=s+0.35_dp*REAL(m-1,dp)
01917     END IF
01918     ! Rule c) if (s,p) shell 0.85 from each electron with n-1, and 1.0
01919     !      from all electrons further in
01920     IF(l1+l2 == 3) THEN
01921        IF(n > 1) THEN
01922           m1=ne(1,n-1)+ne(2,n-1)+ne(3,n-1)+ne(4,n-1)
01923           m2=0
01924           DO i=1,n-2
01925              m2=m2+ne(1,i)+ne(2,i)+ne(3,i)+ne(4,I)
01926           END DO
01927           s=s+0.85_dp*REAL(m1,dp)+1._dp*REAL(m2,dp)
01928        END IF
01929     ELSE
01930        ! Rule d) if (d,f) shell 1.0 from each electron inside
01931        m=0
01932        DO i=1,n-1
01933           m=m+ne(1,i)+ne(2,i)+ne(3,i)+ne(4,i)
01934        END DO
01935        s=s+1._dp*REAL(m,dp)
01936     END IF
01937     ! Slater exponent is (Z-S)/NS
01938     srules = (REAL(z,dp) - s)/xns(n)
01939   END FUNCTION srules
01940 
01941 ! *****************************************************************************
01942 
01943   SUBROUTINE allocate_geminal_basis_set(geminal_basis_set, error)
01944 
01945     ! Allocate a Gaussian-correlated Geminal orbital basis set data set.
01946 
01947     TYPE(geminal_basis_set_type), POINTER    :: geminal_basis_set
01948     TYPE(cp_error_type), INTENT(inout)       :: error
01949 
01950     CHARACTER(len=*), PARAMETER :: routineN = 'allocate_geminal_basis_set', 
01951       routineP = moduleN//':'//routineN
01952 
01953     INTEGER                                  :: stat
01954     LOGICAL                                  :: failure
01955 
01956     failure = .FALSE.
01957     CALL deallocate_geminal_basis_set(geminal_basis_set,error)
01958 
01959     ALLOCATE (geminal_basis_set,STAT=stat)
01960     CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01961 
01962     NULLIFY (geminal_basis_set%cgf_symbol)
01963     NULLIFY (geminal_basis_set%set_radius)
01964     NULLIFY (geminal_basis_set%pgf_radius)
01965     NULLIFY (geminal_basis_set%lmax)
01966     NULLIFY (geminal_basis_set%lmin)
01967     NULLIFY (geminal_basis_set%ls)
01968     NULLIFY (geminal_basis_set%npgf)
01969     NULLIFY (geminal_basis_set%ngem_set)
01970     NULLIFY (geminal_basis_set%nshell)
01971     NULLIFY (geminal_basis_set%l)
01972     NULLIFY (geminal_basis_set%zet)
01973     NULLIFY (geminal_basis_set%zeth)
01974     NULLIFY (geminal_basis_set%first_cgf)
01975     NULLIFY (geminal_basis_set%last_cgf)
01976     NULLIFY (geminal_basis_set%gcc)
01977 
01978   END SUBROUTINE allocate_geminal_basis_set
01979 
01980 ! *****************************************************************************
01981   SUBROUTINE deallocate_geminal_basis_set(geminal_basis_set, error)
01982 
01983     ! Deallocate a Geminal basis set data set.
01984 
01985     TYPE(geminal_basis_set_type), POINTER    :: geminal_basis_set
01986     TYPE(cp_error_type), INTENT(inout)       :: error
01987 
01988     CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_geminal_basis_set', 
01989       routineP = moduleN//':'//routineN
01990 
01991     INTEGER                                  :: stat
01992     LOGICAL                                  :: failure
01993 
01994     failure = .FALSE.
01995     IF (ASSOCIATED(geminal_basis_set))  THEN
01996        IF (ASSOCIATED(geminal_basis_set%cgf_symbol)) THEN
01997           DEALLOCATE (geminal_basis_set%cgf_symbol,STAT=stat)
01998           CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
01999        ENDIF
02000        DEALLOCATE (geminal_basis_set%set_radius,STAT=stat)
02001        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02002        DEALLOCATE (geminal_basis_set%pgf_radius,STAT=stat)
02003        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02004        DEALLOCATE (geminal_basis_set%lmax,STAT=stat)
02005        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02006        DEALLOCATE (geminal_basis_set%lmin,STAT=stat)
02007        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02008        DEALLOCATE (geminal_basis_set%ls,STAT=stat)
02009        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02010        DEALLOCATE (geminal_basis_set%npgf,STAT=stat)
02011        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02012        DEALLOCATE (geminal_basis_set%ngem_set,STAT=stat)
02013        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02014        DEALLOCATE (geminal_basis_set%nshell,STAT=stat)
02015        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02016        DEALLOCATE (geminal_basis_set%l,STAT=stat)
02017        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02018        DEALLOCATE (geminal_basis_set%zet,STAT=stat)
02019        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02020        DEALLOCATE (geminal_basis_set%zeth,STAT=stat)
02021        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02022        DEALLOCATE (geminal_basis_set%gcc,STAT=stat)
02023        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02024        DEALLOCATE (geminal_basis_set%first_cgf,STAT=stat)
02025        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02026        DEALLOCATE (geminal_basis_set%last_cgf,STAT=stat)
02027        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02028        ! now get rid of the full type
02029        DEALLOCATE (geminal_basis_set,STAT=stat)
02030        CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02031     END IF
02032   END SUBROUTINE deallocate_geminal_basis_set
02033 
02034 ! *****************************************************************************
02035   SUBROUTINE get_geminal_basis_set(geminal_basis_set,name,type_restriction,&
02036     cgf_symbol,ngeminals,nset,lmax,lmin,ls,npgf,nshell,l,first_cgf,last_cgf,&
02037     kind_radius,set_radius,pgf_radius,zet,zeth,gcc)
02038 
02039     TYPE(geminal_basis_set_type), POINTER    :: geminal_basis_set
02040     CHARACTER(LEN=default_string_length), 
02041       INTENT(OUT), OPTIONAL                  :: name
02042     CHARACTER(LEN=2), OPTIONAL               :: type_restriction
02043     CHARACTER(LEN=12), DIMENSION(:), 
02044       OPTIONAL, POINTER                      :: cgf_symbol
02045     INTEGER, INTENT(OUT), OPTIONAL           :: ngeminals, nset
02046     INTEGER, DIMENSION(:), OPTIONAL, POINTER :: lmax, lmin, ls, npgf, nshell
02047     INTEGER, DIMENSION(:, :), OPTIONAL, 
02048       POINTER                                :: l, first_cgf, last_cgf
02049     REAL(KIND=dp), OPTIONAL                  :: kind_radius
02050     REAL(KIND=dp), DIMENSION(:), OPTIONAL, 
02051       POINTER                                :: set_radius
02052     REAL(KIND=dp), DIMENSION(:, :), 
02053       OPTIONAL, POINTER                      :: pgf_radius
02054     REAL(KIND=dp), DIMENSION(:, :, :, :), 
02055       OPTIONAL, POINTER                      :: zet, zeth
02056     REAL(KIND=dp), DIMENSION(:, :, :), 
02057       OPTIONAL, POINTER                      :: gcc
02058 
02059     CHARACTER(len=*), PARAMETER :: routineN = 'get_geminal_basis_set', 
02060       routineP = moduleN//':'//routineN
02061 
02062 ! -------------------------------------------------------------------------
02063 
02064     IF (ASSOCIATED(geminal_basis_set)) THEN
02065        IF (PRESENT(name)) name = geminal_basis_set%name
02066        IF (PRESENT(type_restriction)) type_restriction = geminal_basis_set%type_restriction
02067        IF (PRESENT(cgf_symbol)) cgf_symbol => geminal_basis_set%cgf_symbol
02068 
02069        IF (PRESENT(ngeminals)) ngeminals = geminal_basis_set%ngeminals
02070        IF (PRESENT(nset)) nset = geminal_basis_set%nset
02071        IF (PRESENT(npgf)) npgf => geminal_basis_set%npgf
02072        IF (PRESENT(nshell)) nshell => geminal_basis_set%nshell
02073 
02074        IF (PRESENT(lmax)) lmax => geminal_basis_set%lmax
02075        IF (PRESENT(lmin)) lmin => geminal_basis_set%lmin
02076        IF (PRESENT(ls)) ls => geminal_basis_set%ls
02077 
02078        IF (PRESENT(first_cgf)) first_cgf => geminal_basis_set%first_cgf
02079        IF (PRESENT(last_cgf)) last_cgf => geminal_basis_set%last_cgf
02080        IF (PRESENT(l)) l => geminal_basis_set%l
02081 
02082        IF (PRESENT(kind_radius)) kind_radius = geminal_basis_set%kind_radius
02083        IF (PRESENT(set_radius)) set_radius => geminal_basis_set%set_radius
02084        IF (PRESENT(pgf_radius)) pgf_radius => geminal_basis_set%pgf_radius
02085 
02086        IF (PRESENT(zet)) zet => geminal_basis_set%zet
02087        IF (PRESENT(zeth)) zeth => geminal_basis_set%zeth
02088        IF (PRESENT(gcc)) gcc => geminal_basis_set%gcc
02089     ELSE
02090        CALL stop_program(routineN,moduleN,__LINE__,&
02091             "The pointer geminal_basis_set is not associated")
02092     END IF
02093 
02094   END SUBROUTINE get_geminal_basis_set
02095 
02096 ! *****************************************************************************
02097   SUBROUTINE read_geminal_basis_set(element_symbol,basis_set_name,geminal_basis_set,&
02098        para_env,dft_section,basis_section,error)
02099 
02100     CHARACTER(LEN=*), INTENT(IN)             :: element_symbol, basis_set_name
02101     TYPE(geminal_basis_set_type), POINTER    :: geminal_basis_set
02102     TYPE(cp_para_env_type), POINTER          :: para_env
02103     TYPE(section_vals_type), POINTER         :: dft_section
02104     TYPE(section_vals_type), OPTIONAL, 
02105       POINTER                                :: basis_section
02106     TYPE(cp_error_type), INTENT(inout)       :: error
02107 
02108     CHARACTER(len=*), PARAMETER :: routineN = 'read_geminal_basis_set', 
02109       routineP = moduleN//':'//routineN
02110 
02111     CHARACTER(len=20*default_string_length)  :: line_att
02112     CHARACTER(LEN=240)                       :: line
02113     CHARACTER(LEN=242)                       :: line2
02114     CHARACTER(len=default_path_length)       :: basis_set_file_name
02115     CHARACTER(LEN=LEN(basis_set_name))       :: bsname
02116     CHARACTER(LEN=LEN(basis_set_name)+2)     :: bsname2
02117     CHARACTER(LEN=LEN(element_symbol))       :: symbol
02118     CHARACTER(LEN=LEN(element_symbol)+2)     :: symbol2
02119     INTEGER :: i, ico, ipgf, irep, iset, ishell, istr, lshell, lss, lval, 
02120       maxco, maxl, maxpgf, maxshell, ncgf, ngeminals, nset, stat, strlen1, 
02121       strlen2
02122     INTEGER, DIMENSION(:), POINTER           :: lmax, lmin, ls, npgf, nshell
02123     INTEGER, DIMENSION(:, :), POINTER        :: l
02124     LOGICAL                                  :: failure, found, is_ok, isrr, 
02125                                                 isrs, match, read_from_input
02126     REAL(KIND=dp)                            :: aa, bb, cc
02127     REAL(KIND=dp), DIMENSION(2, 2)           :: zl
02128     REAL(KIND=dp), DIMENSION(:, :, :), 
02129       POINTER                                :: gcc
02130     REAL(KIND=dp), DIMENSION(:, :, :, :), 
02131       POINTER                                :: zet
02132     TYPE(cp_parser_type), POINTER            :: parser
02133     TYPE(cp_sll_val_type), POINTER           :: list
02134     TYPE(val_type), POINTER                  :: val
02135 
02136     failure = .FALSE.
02137     line = ""
02138     line2 = ""
02139     symbol = ""
02140     symbol2 = ""
02141     bsname = ""
02142     bsname2 = ""
02143 
02144     IF ( basis_set_name == "" ) THEN
02145       geminal_basis_set%name = element_symbol//"_Geminal"
02146     ELSE
02147       geminal_basis_set%name = basis_set_name
02148     END IF
02149     read_from_input = .FALSE.
02150     IF (PRESENT(basis_section)) THEN
02151        CALL section_vals_get(basis_section,explicit=read_from_input, error=error)
02152     END IF
02153     IF (.NOT.read_from_input) THEN
02154        CALL section_vals_val_get(dft_section,"GEMINAL_FILE_NAME",&
02155             c_val=basis_set_file_name,error=error)
02156        NULLIFY(parser)
02157        CALL parser_create(parser,basis_set_file_name,para_env=para_env,error=error)
02158     END IF
02159 
02160     ! Search for the requested basis set in the basis set file
02161     ! until the basis set is found or the end of file is reached
02162 
02163     bsname = basis_set_name
02164     symbol = element_symbol
02165     irep   = 0
02166     match = .FALSE.
02167 
02168     nset = 0
02169     maxshell=0
02170     maxpgf=0
02171     maxco=0
02172     ngeminals=0
02173     geminal_basis_set%nset = nset
02174     geminal_basis_set%ngeminals = ngeminals
02175     CALL reallocate(geminal_basis_set%lmax,1,nset)
02176     CALL reallocate(geminal_basis_set%lmin,1,nset)
02177     CALL reallocate(geminal_basis_set%ls,1,nset)
02178     CALL reallocate(geminal_basis_set%nshell,1,nset)
02179     CALL reallocate(geminal_basis_set%l,1,maxshell,1,nset)
02180     CALL reallocate(geminal_basis_set%zet,1,2,1,2,1,maxpgf,1,nset)
02181     CALL reallocate(geminal_basis_set%gcc,1,maxpgf,1,maxshell,1,nset)
02182     CALL reallocate(geminal_basis_set%set_radius,1,nset)
02183     CALL reallocate(geminal_basis_set%pgf_radius,1,maxpgf,1,nset)
02184     CALL reallocate(geminal_basis_set%first_cgf,1,maxshell,1,nset)
02185     CALL reallocate(geminal_basis_set%last_cgf,1,maxshell,1,nset)
02186 
02187     search_loop: DO
02188 
02189          IF (read_from_input) THEN
02190             NULLIFY(list,val)
02191             found = .TRUE.
02192             CALL section_vals_list_get(basis_section,"_DEFAULT_KEYWORD_",list=list,error=error)
02193          ELSE
02194             CALL parser_search_string(parser,TRIM(bsname),.TRUE.,found,line,error=error)
02195          END IF
02196          IF (found) THEN
02197             CALL uppercase(symbol)
02198             CALL uppercase(bsname)
02199 
02200             IF (read_from_input) THEN
02201                match = .TRUE.
02202             ELSE
02203                match = .FALSE.
02204                CALL uppercase(line)
02205                ! Check both the element symbol and the basis set name
02206                line2 = " "//line//" "
02207                symbol2 = " "//TRIM(symbol)//" "
02208                bsname2 = " "//TRIM(bsname)//" "
02209                strlen1 = LEN_TRIM(symbol2) + 1
02210                strlen2 = LEN_TRIM(bsname2) + 1
02211 
02212                IF ( (INDEX(line2,symbol2(:strlen1)) > 0).AND.&
02213                     (INDEX(line2,bsname2(:strlen2)) > 0) ) match = .TRUE.
02214             END IF
02215             IF (match) THEN
02216                NULLIFY (gcc,l,lmax,lmin,npgf,nshell,zet,ls)
02217                ! Read the basis set information
02218                IF (read_from_input) THEN
02219                   is_ok=cp_sll_val_next(list,val,error=error)
02220                   IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,&
02221                        __LINE__,&
02222                        "Error reading the Geminal basis set from input file!!")
02223                   CALL val_get(val,c_val=line_att,error=error)
02224                   READ(line_att,*)nset
02225                ELSE
02226                   CALL parser_get_object(parser,nset,newline=.TRUE.,error=error)
02227                   IF (PRESENT(basis_section)) THEN
02228                      irep = irep + 1
02229                      WRITE(line_att,'(1X,I0)')nset
02230                      CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
02231                           c_val=TRIM(line_att), error=error)
02232                   END IF
02233                END IF
02234 
02235                CALL reallocate(npgf,1,nset)
02236                CALL reallocate(nshell,1,nset)
02237                CALL reallocate(lmax,1,nset)
02238                CALL reallocate(lmin,1,nset)
02239                CALL reallocate(ls,1,nset)
02240 
02241                maxl = 0
02242                maxpgf = 0
02243                maxshell = 0
02244 
02245                DO iset=1,nset
02246                   IF (read_from_input) THEN
02247                      is_ok=cp_sll_val_next(list,val,error=error)
02248                      IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,&
02249                           __LINE__,&
02250                           "Error reading the Geminal basis set from input file!!")
02251                      CALL val_get(val,c_val=line_att,error=error)
02252                      READ(line_att,*)lmin(iset)
02253                      CALL remove_word(line_att)
02254                      READ(line_att,*)lmax(iset)
02255                      CALL remove_word(line_att)
02256                      READ(line_att,*)ls(iset)
02257                      CALL remove_word(line_att)
02258                      READ(line_att,*)npgf(iset)
02259                      CALL remove_word(line_att)
02260                   ELSE
02261                      line_att = ""
02262                      CALL parser_get_object(parser,lmin(iset),newline=.TRUE.,error=error)
02263                      CALL parser_get_object(parser,lmax(iset),error=error)
02264                      CALL parser_get_object(parser,ls(iset),error=error)
02265                      CALL parser_get_object(parser,npgf(iset),error=error)
02266                      WRITE(line_att,'(4(1X,I0))')lmin(iset),lmax(iset),ls(iset),npgf(iset)
02267                   END IF
02268                   maxl = MAX(maxl,lmax(iset),ls(iset))
02269                   IF (npgf(iset) > maxpgf) THEN
02270                      maxpgf = npgf(iset)
02271                      CALL reallocate(zet,1,2,1,2,1,maxpgf,1,nset)
02272                      CALL reallocate(gcc,1,maxpgf,1,maxshell,1,nset)
02273                   END IF
02274                   nshell(iset) = 0
02275                   DO lshell=lmin(iset),lmax(iset)
02276                      IF (read_from_input) THEN
02277                         READ(line_att,*)ishell
02278                         CALL remove_word(line_att)
02279                      ELSE
02280                         CALL parser_get_object(parser,ishell,error=error)
02281                         istr=LEN_TRIM(line_att)+1
02282                         WRITE(line_att(istr:),'(1X,I0)')ishell
02283                      END IF
02284                      nshell(iset) = nshell(iset) + ishell
02285                      IF (nshell(iset) > maxshell) THEN
02286                         maxshell = nshell(iset)
02287                         CALL reallocate(l,1,maxshell,1,nset)
02288                         CALL reallocate(gcc,1,maxpgf,1,maxshell,1,nset)
02289                      END IF
02290                      DO i=1,ishell
02291                         l(nshell(iset)-ishell+i,iset) = lshell
02292                      END DO
02293                   END DO
02294                   IF (.NOT.read_from_input) THEN
02295                      IF (PRESENT(basis_section)) THEN
02296                         irep = irep + 1
02297                         CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
02298                              c_val=TRIM(line_att), error=error)
02299                         line_att = ""
02300                      END IF
02301                   ELSE
02302                      IF (LEN_TRIM(line_att)/=0) CALL stop_program(routineN,moduleN,&
02303                           __LINE__,&
02304                           "Error reading the Geminal basis from input file!!")
02305                   END IF
02306                   DO ipgf=1,npgf(iset)
02307                      IF (read_from_input) THEN
02308                         is_ok=cp_sll_val_next(list,val,error=error)
02309                         IF (.NOT.is_ok) CALL stop_program(routineN,moduleN,&
02310                              __LINE__,&
02311                              "Error reading the Basis set from input file!!")
02312                         CALL val_get(val,c_val=line_att,error=error)
02313                         READ(line_att,*)zl(1,1),zl(2,2),zl(1,2),(gcc(ipgf,ishell,iset),ishell=1,nshell(iset))
02314                         zet(1,1,ipgf,iset)=zl(1,1)**2
02315                         zet(2,2,ipgf,iset)=zl(2,2)**2 + zl(1,2)**2
02316                         zet(1,2,ipgf,iset)=zl(1,1)*zl(1,2)
02317                         zet(2,1,ipgf,iset)=zl(1,1)*zl(1,2)
02318                      ELSE
02319                         CALL parser_get_object(parser,zl(1,1),newline=.TRUE.,error=error)
02320                         CALL parser_get_object(parser,zl(2,2),error=error)
02321                         CALL parser_get_object(parser,zl(1,2),error=error)
02322                         zet(1,1,ipgf,iset)=zl(1,1)**2
02323                         zet(2,2,ipgf,iset)=zl(2,2)**2 + zl(1,2)**2
02324                         zet(1,2,ipgf,iset)=zl(1,1)*zl(1,2)
02325                         zet(2,1,ipgf,iset)=zl(1,1)*zl(1,2)
02326                         DO ishell=1,nshell(iset)
02327                            CALL parser_get_object(parser,gcc(ipgf,ishell,iset),error=error)
02328                         END DO
02329                         IF (PRESENT(basis_section)) THEN
02330                            irep = irep + 1
02331                            WRITE(line_att,'(100E24.16)')zet(1,1,ipgf,iset),zet(2,2,ipgf,iset),zet(1,2,ipgf,iset),&
02332                                 (gcc(ipgf,ishell,iset),ishell=1,nshell(iset))
02333                            CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
02334                                 c_val=TRIM(line_att), error=error)
02335                         END IF
02336                      END IF
02337                   END DO
02338                END DO
02339 
02340                ! Maximum angular momentum quantum number of the atomic kind
02341 
02342                CALL init_orbital_pointers(2*maxl)
02343 
02344                ! Allocate the global variables
02345 
02346                geminal_basis_set%nset = nset
02347                CALL reallocate(geminal_basis_set%lmax,1,nset)
02348                CALL reallocate(geminal_basis_set%lmin,1,nset)
02349                CALL reallocate(geminal_basis_set%ls,1,nset)
02350                CALL reallocate(geminal_basis_set%nshell,1,nset)
02351                CALL reallocate(geminal_basis_set%npgf,1,nset)
02352                CALL reallocate(geminal_basis_set%ngem_set,1,nset)
02353                CALL reallocate(geminal_basis_set%l,1,maxshell,1,nset)
02354                CALL reallocate(geminal_basis_set%zet,1,2,1,2,1,maxpgf,1,nset)
02355                CALL reallocate(geminal_basis_set%zeth,1,2,1,2,1,maxpgf,1,nset)
02356                CALL reallocate(geminal_basis_set%gcc,1,maxpgf,1,maxshell,1,nset)
02357 
02358                ! Copy the basis set information into the data structure
02359 
02360                DO iset=1,nset
02361                   geminal_basis_set%lmax(iset) = lmax(iset)
02362                   geminal_basis_set%lmin(iset) = lmin(iset)
02363                   geminal_basis_set%ls(iset)   = ls(iset)
02364                   geminal_basis_set%npgf(iset) = npgf(iset)
02365                   geminal_basis_set%nshell(iset) = nshell(iset)
02366                   DO ishell=1,nshell(iset)
02367                      geminal_basis_set%l(ishell,iset) = l(ishell,iset)
02368                      DO ipgf=1,npgf(iset)
02369                         geminal_basis_set%gcc(ipgf,ishell,iset) = gcc(ipgf,ishell,iset)
02370                      END DO
02371                   END DO
02372                   DO ipgf=1,npgf(iset)
02373                      geminal_basis_set%zet(1:2,1:2,ipgf,iset) = zet(1:2,1:2,ipgf,iset)
02374                      geminal_basis_set%zeth(1,1,ipgf,iset) = zet(1,1,ipgf,iset)
02375                      geminal_basis_set%zeth(2,2,ipgf,iset) = zet(2,2,ipgf,iset)
02376                      geminal_basis_set%zeth(1,2,ipgf,iset) = -zet(1,2,ipgf,iset)
02377                      geminal_basis_set%zeth(2,1,ipgf,iset) = -zet(2,1,ipgf,iset)
02378                   END DO
02379                END DO
02380 
02381                ! Initialise the depending atomic kind information
02382 
02383                CALL reallocate(geminal_basis_set%set_radius,1,nset)
02384                CALL reallocate(geminal_basis_set%pgf_radius,1,maxpgf,1,nset)
02385                CALL reallocate(geminal_basis_set%first_cgf,1,maxshell,1,nset)
02386                CALL reallocate(geminal_basis_set%last_cgf,1,maxshell,1,nset)
02387                CALL reallocate(geminal_basis_set%ngem_set,1,nset)
02388 
02389                maxco = 0
02390                ncgf = 0
02391 
02392                DO iset=1,nset
02393                   geminal_basis_set%ngem_set(iset) = 0
02394                   DO ishell=1,nshell(iset)
02395                      lshell = geminal_basis_set%l(ishell,iset)
02396                      geminal_basis_set%first_cgf(ishell,iset) = ncgf + 1
02397                      ncgf = ncgf + nco(lshell)*nco(ls(iset))
02398                      geminal_basis_set%last_cgf(ishell,iset) = ncgf
02399                      geminal_basis_set%ngem_set(iset) =&
02400                           geminal_basis_set%ngem_set(iset) + nco(lshell)*nco(ls(iset))
02401                   END DO
02402                   maxco = MAX(maxco,npgf(iset)*ncoset(lmax(iset))*ncoset(ls(iset)))
02403                END DO
02404 
02405                geminal_basis_set%ngeminals = ncgf
02406 
02407                ALLOCATE (geminal_basis_set%cgf_symbol(ncgf),STAT=stat)
02408                CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02409 
02410                ncgf = 0
02411 
02412                DO iset=1,nset
02413                   lval = geminal_basis_set%ls(iset)
02414                   DO ishell=1,nshell(iset)
02415                      lshell = geminal_basis_set%l(ishell,iset)
02416                      DO ico=ncoset(lshell-1)+1,ncoset(lshell)
02417                         DO lss=1,nco(lval)
02418                            ncgf = ncgf + 1
02419                            geminal_basis_set%cgf_symbol(ncgf) =&
02420                                 cgf_symbol(lshell+1,indco(1:3,ico))//"*"//cgf_symbol(lval+1,indco(1:3,lss))
02421                         END DO
02422                      END DO
02423                   END DO
02424                END DO
02425 
02426                DEALLOCATE (gcc,lmax,lmin,ls,l,npgf,nshell,zet,STAT=stat)
02427                CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
02428 
02429                EXIT search_loop
02430             END IF
02431          ELSE
02432             ! Stop program, if the end of file is reached
02433             CALL stop_program(routineN,moduleN,__LINE__,&
02434                  "The requested geminal basis set <"//TRIM(bsname)//&
02435                  "> for element <"//TRIM(symbol)//"> was not "//&
02436                  "found in the basis set file <"//&
02437                  TRIM(basis_set_file_name)//">",para_env)
02438          END IF
02439     END DO search_loop
02440 
02441     ! is this a RS or RR basis
02442     isrs = .TRUE.
02443     isrr = .TRUE.
02444     DO iset=1,geminal_basis_set%nset
02445        DO ipgf=1,geminal_basis_set%npgf(iset)
02446          aa = geminal_basis_set%zet(1,1,ipgf,iset)
02447          bb = geminal_basis_set%zet(2,2,ipgf,iset)
02448          cc = geminal_basis_set%zet(1,2,ipgf,iset)
02449          IF (ABS(cc) > 1.E-6) isrs = .FALSE.
02450          IF (ABS(aa-bb) > 1.E-6 .OR. aa < ABS(cc) ) isrr = .FALSE.
02451        END DO
02452     END DO
02453     IF ( isrs ) THEN
02454       geminal_basis_set%type_restriction = "RS"
02455     ELSEIF ( isrr ) THEN
02456       geminal_basis_set%type_restriction = "RR"
02457     ELSE
02458       geminal_basis_set%type_restriction = "NO"
02459     END IF
02460 
02461     IF (.NOT.read_from_input) THEN
02462        CALL parser_release(parser,error=error)
02463        IF ((match).AND.(PRESENT(basis_section))) THEN
02464           ! Dump the read basis set in the basis section
02465           irep = irep + 1
02466           WRITE(line_att,'(A)')"         # Geminal basis set name: "//bsname2(:strlen2)//&
02467                &" for symbol: "//symbol2(:strlen1)
02468           CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
02469                c_val=TRIM(line_att), error=error)
02470           irep = irep + 1
02471           WRITE(line_att,'(A)')"         # Geminal basis set read from the basis set filename: "//&
02472                &TRIM(basis_set_file_name)
02473           CALL section_vals_val_set(basis_section,"_DEFAULT_KEYWORD_",i_rep_val=irep,&
02474                c_val=TRIM(line_att), error=error)
02475        END IF
02476     END IF
02477 
02478   END SUBROUTINE read_geminal_basis_set
02479 
02480 ! *****************************************************************************
02481   SUBROUTINE set_geminal_basis_set(geminal_basis_set,name,type_restriction,&
02482                                    kind_radius,set_radius,pgf_radius)
02483 
02484     TYPE(geminal_basis_set_type), POINTER    :: geminal_basis_set
02485     CHARACTER(LEN=default_string_length), 
02486       INTENT(IN), OPTIONAL                   :: name
02487     CHARACTER(LEN=2), OPTIONAL               :: type_restriction
02488     REAL(KIND=dp), OPTIONAL                  :: kind_radius
02489     REAL(KIND=dp), DIMENSION(:), OPTIONAL, 
02490       POINTER                                :: set_radius
02491     REAL(KIND=dp), DIMENSION(:, :), 
02492       OPTIONAL, POINTER                      :: pgf_radius
02493 
02494     CHARACTER(len=*), PARAMETER :: routineN = 'set_geminal_basis_set', 
02495       routineP = moduleN//':'//routineN
02496 
02497 ! -------------------------------------------------------------------------
02498 
02499     IF (ASSOCIATED(geminal_basis_set)) THEN
02500        IF (PRESENT(name)) geminal_basis_set%name = name
02501        IF (PRESENT(type_restriction)) geminal_basis_set%type_restriction = type_restriction
02502 
02503        IF (PRESENT(kind_radius)) geminal_basis_set%kind_radius = kind_radius
02504        IF (PRESENT(set_radius)) geminal_basis_set%set_radius(:) = set_radius(:)
02505        IF (PRESENT(pgf_radius)) geminal_basis_set%pgf_radius(:,:) = pgf_radius(:,:)
02506 
02507     ELSE
02508        CALL stop_program(routineN,moduleN,__LINE__,&
02509             "The pointer geminal_basis_set is not associated")
02510     END IF
02511 
02512   END SUBROUTINE set_geminal_basis_set
02513 
02514 ! *****************************************************************************
02515   SUBROUTINE write_geminal_basis_set(geminal_basis_set,output_unit,error)
02516 
02517     TYPE(geminal_basis_set_type), POINTER    :: geminal_basis_set
02518     INTEGER, INTENT(in)                      :: output_unit
02519     TYPE(cp_error_type), INTENT(inout)       :: error
02520 
02521     CHARACTER(len=*), PARAMETER :: routineN = 'write_geminal_basis_set', 
02522       routineP = moduleN//':'//routineN
02523 
02524     INTEGER                                  :: ipgf, iset, ishell
02525 
02526 ! -------------------------------------------------------------------------
02527 
02528     IF (ASSOCIATED(geminal_basis_set).AND.(output_unit>0)) THEN
02529        WRITE (UNIT=output_unit,FMT="(/,T3,A,T51,A30)")&
02530             "Gaussian geminal basis set information for",TRIM(geminal_basis_set%name)
02531        WRITE (UNIT=output_unit,FMT="(T5,A,T79,A2)")&
02532             "Geminal type restrictions",TRIM(geminal_basis_set%type_restriction)
02533        WRITE (UNIT=output_unit,FMT="(T5,A,T45,A)")&
02534             "Set  Shell   l(R) l(S)","Exponents                Coefficient"
02535 
02536        DO iset=1,geminal_basis_set%nset
02537           DO ishell=1,geminal_basis_set%nshell(iset)
02538              WRITE (UNIT=output_unit,&
02539                   FMT="(T5,I3,2X,I4,I5,I5,(T29,4F13.6))")&
02540                   iset,ishell,&
02541                   geminal_basis_set%l(ishell,iset),&
02542                   geminal_basis_set%ls(iset),&
02543                   (geminal_basis_set%zet(1,1,ipgf,iset),&
02544                   geminal_basis_set%zet(2,2,ipgf,iset),&
02545                   geminal_basis_set%zet(1,2,ipgf,iset),&
02546                   geminal_basis_set%gcc(ipgf,ishell,iset),&
02547                   ipgf=1,geminal_basis_set%npgf(iset))
02548           END DO
02549        END DO
02550 
02551     END IF
02552 
02553   END SUBROUTINE write_geminal_basis_set
02554 ! *****************************************************************************
02555 
02556 END MODULE basis_set_types