diff --git a/lioamber/Makefile.depends b/lioamber/Makefile.depends index 717222bac..f30c81c2d 100644 --- a/lioamber/Makefile.depends +++ b/lioamber/Makefile.depends @@ -180,13 +180,9 @@ $(tmplist:%.o=$(OBJPATH)/%.o) : $(OBJPATH)/debug_tools.mod ################################################################################ OBJECTS += basis_data.o -OBJECTS += basis_subs.o tmplist := -tmplist += ehrensubs.o init_lio.o -$(tmplist:%.o=$(OBJPATH)/%.o) : $(OBJPATH)/basis_subs.mod - -tmplist := basis_subs.o +tmplist += ehrensubs.o init_lio.o lio_finalize.o $(tmplist:%.o=$(OBJPATH)/%.o) : $(OBJPATH)/basis_data.mod @@ -385,7 +381,7 @@ SRCDIRS += maskrmm OBJECTS += maskrmm.o tmplist := -tmplist += ehrensubs.o basis_subs.o +tmplist += ehrensubs.o $(tmplist:%.o=$(OBJPATH)/%.o) : $(OBJPATH)/maskrmm.mod # diff --git a/lioamber/basis_data.f90 b/lioamber/basis_data.f90 index f99a9444e..63654384f 100644 --- a/lioamber/basis_data.f90 +++ b/lioamber/basis_data.f90 @@ -1,15 +1,1073 @@ !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! +!%% BASIS.F90 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! +! This file contains two modules: basis_data, containing basis set function ! +! data, and basis_subs, containing the following subroutines: ! +! ! +! External access: ! +! · basis_init() : Initializes basis data. ! +! --> Requires basis and fitting sets names (or basis file name), number of ! +! atoms and nuclear atomic charges. ! +! · basis_setup_ehren(): Initialises data for ehrenfest runs ! +! --> No inputs required. ! +! · basis_deinit() : Deallocates basis data. ! +! --> No inputs required. ! +! ! +! Intended for basis_subs internal-only access: ! +! · basis_read_internal(): Reads basis functions from /dat directory. ! +! · basis_read_external(): Reads basis functions form a custom basis file. ! +! · basis_set_size() : Prereads data and sets array sizes ! +! · check_basis() : Checks if all atoms have basis functions assigned. ! +! ! +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! + +! TODO: INVERT C and A matrices, since loops are faster that way. In Ehrenfest +! these are already inverted. + module basis_data implicit none - integer :: basis_size - integer :: basis_size_s - integer :: basis_size_p - integer :: basis_size_d - integer :: maximum_contractions - integer, allocatable, dimension(:) :: orbital_contractions - integer, allocatable, dimension(:) :: parent_atom - integer, allocatable, dimension(:,:) :: angular_momentum - real*8 , allocatable, dimension(:,:) :: gauss_expo - real*8 , allocatable, dimension(:,:) :: gauss_coef -end module + ! Namelist inputfiles + ! int_basis : Use LIO internal basis (located in /dat) + ! rMax : Maximum exponent for double-precision integrals. + ! rMaxs : Maximum exponent for single-precision integrals. + ! norm : Normalize integrals (deprecated). + logical :: int_basis = .true. + double precision :: rMax = 16.0D0 + double precision :: rMaxs = 5.0D0 + logical :: norm = .true. + + ! Single variables + ! M : number of basis functions. + ! Md : number of fitting functions. + ! kknums: number of single-precision two-center integrals. + ! kknumd: number of double-precision two-center integrals. + integer :: M = 0 + integer :: Md = 0 + integer :: kknums = 0 + integer :: kknumd = 0 + + ! Preset arrays + ! nShell : Number of basis functions for each shell (s,p,d,f) + ! nShelld: Number of auxiliary basis functions for each shell (s,p,d,f) + ! max_c_per_atom: Maximum number of contractions for a single function. + ! max_f_per_atom: Maximum number of functions for a single atom. + integer :: nShell(0:3) + integer :: nShelld(0:3) + integer :: max_c_per_atom + integer :: max_f_per_atom + + ! Allocatable arrays + ! Nuc(i) : Index of the atom containing function i. + ! Nucd(i) : Index of the atom containing fitting function i. + ! nCont(i) : Number of contractions for function i. + ! nContd(i) : Number of contractions for fitting function i. + ! ang_mom(i) : Angular momentum of function i. Or an angry mother. Your choice. + ! ang_momd(i): Angular momentum of auxiliary function i. + ! a(i,j) : Exponent for function i, contraction j. + ! c(i,j) : Coefficient for function i, contraction j. + ! ad(i,j) : Exponent for auxiliary function i, contraction j. + ! cd(i,j) : Coefficient for auxiliary function i, contraction j. + ! atmin(i) : The minimum exponent found for atom i. + ! kkInd(:) : Index for double-precision two-center integrals. + ! kkInds(:) : Index for single-precision two-center integrals. + ! cool(:) : Storage for two-center integrals in double precision. + ! cools(:) : Storage for two-center integrals in single precision. + integer , allocatable :: Nuc(:) + integer , allocatable :: Nucd(:) + integer , allocatable :: nCont(:) + integer , allocatable :: nContd(:) + integer , allocatable :: ang_mom(:) + integer , allocatable :: ang_momd(:) + integer , allocatable :: kkInd(:) + integer , allocatable :: kkInds(:) + integer , allocatable :: natomc(:) + integer , allocatable :: nns(:) + integer , allocatable :: nnp(:) + integer , allocatable :: nnd(:) + integer , allocatable :: nnps(:) + integer , allocatable :: nnpd(:) + integer , allocatable :: nnpp(:) + integer , allocatable :: jatc(:,:) + double precision, allocatable :: a(:,:) + double precision, allocatable :: c(:,:) + double precision, allocatable :: ad(:,:) + double precision, allocatable :: cd(:,:) + double precision, allocatable :: af(:) + double precision, allocatable :: atmin(:) + double precision, allocatable :: cool(:) + real , allocatable :: cools(:) + + + ! Temporary for EHRENFEST + double precision, allocatable :: a_ehren(:,:) + double precision, allocatable :: c_ehren(:,:) + integer , allocatable :: ang_mom_ehren(:,:) + + ! GLOBAL PARAMETERS + ! Degeneracy for each angular momentum + integer , parameter :: ANG_DEG(0:3) = (/1, 3, 6, 10/) + integer , parameter :: MAX_CONTRACT = 50 +contains +end module basis_data +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! +module basis_subs + implicit none + integer, parameter :: TMP_OPEN_UID = 9999 +contains + +subroutine basis_init(basis_name, fitting_name, n_atoms, atom_Z, out_stat) + !use basis_data, only: M, Md, int_basis, Nuc, Nucd, nCont, nContd, a, c, ad, & + use garcha_mod, only: M, Md, int_basis, Nuc, Nucd, nCont, nContd, a, c, ad, & + cd, atmin, nns, nnp, nnd, nshell, nshelld, norm, af, & + indexii, indexiid, natomc, jatc, nnps, nnpp, nnpd + use basis_data, only: ang_mom, ang_momd, max_f_per_atom, max_c_per_atom + implicit none + ! Inputs: + ! n_atoms : the number of atoms in the QM system. + ! atom_Z(n_atoms): nuclear atomic charge of said atoms. + ! basis_name : the name of the basis set or basis file. + ! fitting_name : the name of the fitting set (not used when int_basis is + ! true). + ! Outputs: + ! out_stat : Returns an iostat-like value for error handling. + + integer , intent(in) :: n_atoms, atom_Z(n_atoms) + integer , intent(out) :: out_stat + character(len=*), intent(inout) :: basis_name, fitting_name + + integer :: icount, iostat + integer, allocatable :: atom_count(:) + logical, allocatable :: atom_basis_chk(:), atom_fitting_chk(:) + + out_stat = 0 + allocate(atom_count(0:120), atom_basis_chk(0:120), atom_fitting_chk(0:120)) + + ! Precalculates the amount of atoms of a certain type, in order to + ! correctly asses the sizes of M and Md + atom_count = 0 + do icount = 1, n_atoms + atom_count(atom_Z(icount)) = atom_count(atom_Z(icount)) + 1 + enddo + + call basis_set_size(M, Md, max_f_per_atom, max_c_per_atom, atom_basis_chk, & + atom_fitting_chk, basis_name, fitting_name, atom_count, & + atom_Z, n_atoms, int_basis, iostat) + if (iostat .gt. 0) then + out_stat = 1 + return + endif + + call check_basis(atom_basis_chk, atom_fitting_chk, n_atoms, atom_Z, iostat) + if (iostat .gt. 0) then + out_stat = 2 + return + endif + + allocate(c(M, max_c_per_atom), a(M, max_c_per_atom), cd(Md, max_c_per_atom),& + ad(Md, max_c_per_atom), nCont(M), nContd(Md), ang_mom(M), & + ang_momd(Md), Nuc(M), Nucd(Md), indexii(M), indexiid(Md), af(Md)) + allocate(atmin(n_atoms), nns(n_atoms), nnp(n_atoms), nnd(n_atoms), & + nnps(n_atoms), nnpp(n_atoms), nnpd(n_atoms), natomc(n_atoms), & + jatc(n_atoms,n_atoms)) + + ! Initializes everything to 0. + c = 0.0D0 ; a = 0.0D0; nCont = 0; ang_mom = 0; nuc = 0; af = 0.0D0 + cd = 0.0D0 ; ad = 0.0D0; nContd = 0; ang_momd = 0; nucd = 0 + + nns = 0 ; nnp = 0 ; nnd = 0 ; natomc = 0 + nnps = 0 ; nnpp = 0 ; nnpd = 0 ; jatc = 0 + + + if (int_basis) then + call read_basis_internal(basis_name, fitting_name, M, Md, n_atoms, norm, & + max_f_per_atom, max_c_per_atom, atom_Z, c, a, & + cd, ad, nCont, nContd, ang_mom, ang_momd, Nuc, & + Nucd, atmin, nns, nnp, nnd, nShell, nShelld, & + iostat) + else + call read_basis_external(basis_name, M,Md, n_atoms, norm, max_f_per_atom,& + max_c_per_atom, atom_Z, c, a, cd, ad, nCont, & + nContd, ang_mom, ang_momd, Nuc, Nucd, atmin, & + nns, nnp, nnd, nShell, nShelld, iostat) + endif + if (iostat .gt. 0) then + out_stat = 3 + return + endif + + ! Reorders basis: first all s, then all p, then all d. + call reorder_basis(a, c, Nuc, nCont, indexii, M, max_c_per_atom, & + ang_mom, nShell) + call reorder_basis(ad, cd, Nucd, nContd, indexiid, Md, max_c_per_atom, & + ang_momd, nShelld) + + deallocate(atom_count, atom_basis_chk, atom_fitting_chk, ang_mom, ang_momd) +end subroutine basis_init + +subroutine basis_deinit() + !use basis_data, only: Nuc, Nucd, nCont, nContd, a, c, ad, cd, atmin, nns, & + use garcha_mod , only: Nuc, Nucd, nCont, nContd, a, c, ad, cd, atmin, nns, & + nnp, nnd, af, indexii, indexiid, natomc, jatc, nnps,& + nnpp, nnpd + + implicit none + + ! M or Md sized. + if (allocated(c)) deallocate(c) + if (allocated(a)) deallocate(a) + if (allocated(cd)) deallocate(cd) + if (allocated(ad)) deallocate(ad) + if (allocated(af)) deallocate(af) + if (allocated(ncont)) deallocate(ncont) + if (allocated(ncontd)) deallocate(ncontd) + if (allocated(nuc)) deallocate(nuc) + if (allocated(nucd)) deallocate(nucd) + if (allocated(indexii)) deallocate(indexii) + if (allocated(indexiid)) deallocate(indexiid) + + ! natom sized. + if (allocated(atmin)) deallocate(atmin) + if (allocated(jatc)) deallocate(jatc) + if (allocated(natomc)) deallocate(natomc) + if (allocated(nns)) deallocate(nns) + if (allocated(nnp)) deallocate(nnp) + if (allocated(nnd)) deallocate(nnd) + if (allocated(nnps)) deallocate(nnps) + if (allocated(nnpp)) deallocate(nnpp) + if (allocated(nnpd)) deallocate(nnpd) + +end subroutine basis_deinit + +subroutine basis_setup_ehren() + use basis_data, only: a_ehren, c_ehren, ang_mom_ehren, max_c_per_atom + use garcha_mod, only: a, c, nShell, M + + implicit none + integer :: icount + + ! Transposes coefficient and exponent matrices. + allocate (c_ehren(max_c_per_atom,M), a_ehren(max_c_per_atom, M)) + do icount = 1, M + a_ehren(:,icount) = a(icount,:) + c_ehren(:,icount) = c(icount,:) + enddo + + ! Calculates angular momenta and directions + allocate(ang_mom_ehren(3,M)) + ang_mom_ehren = 0 + do icount = nShell(0)+1, nShell(0)+nShell(1), 3 + ang_mom_ehren(1, icount ) = 1 + ang_mom_ehren(2, icount+1) = 1 + ang_mom_ehren(3, icount+2) = 1 + enddo + do icount = nShell(0)+nShell(1)+1, nShell(0)+nShell(1)+nShell(2), 6 + ang_mom_ehren(1,icount+0) = 2 ! dxx (x) + ang_mom_ehren(1,icount+1) = 1 ! dxy (x) + ang_mom_ehren(2,icount+1) = 1 ! dxy (y) + ang_mom_ehren(2,icount+2) = 2 ! dyy (y) + ang_mom_ehren(1,icount+3) = 1 ! dxz (x) + ang_mom_ehren(3,icount+3) = 1 ! dxz (z) + ang_mom_ehren(2,icount+4) = 1 ! dyz (y) + ang_mom_ehren(3,icount+4) = 1 ! dyz (z) + ang_mom_ehren(3,icount+5) = 2 ! dzz (z) + + c_ehren(:,icount ) = c_ehren(:,icount ) / dsqrt(3.0D0) + c_ehren(:,icount+2) = c_ehren(:,icount+2) / dsqrt(3.0D0) + c_ehren(:,icount+5) = c_ehren(:,icount+5) / dsqrt(3.0D0) + enddo +end subroutine basis_setup_ehren + +subroutine basis_set_size(basis_size, aux_size, max_f_per_atom, max_c_per_atom,& + atom_bas_done, atom_fit_done, basis_file, & + fitting_file, atom_count, atom_Z, n_atoms, & + use_internal, iostatus) + use basis_data, only: ANG_DEG, MAX_CONTRACT + implicit none + integer , intent(in) :: n_atoms, atom_Z(n_atoms), & + atom_count(0:120) + logical , intent(in) :: use_internal + integer , intent(out) :: basis_size, aux_size, max_f_per_atom, & + max_c_per_atom, iostatus + logical , intent(out) :: atom_bas_done(0:120), & + atom_fit_done(0:120) + character(len=*), intent(inout) :: basis_file, fitting_file + + logical :: file_exists + integer :: file_uid = TMP_OPEN_UID + integer :: file_iostat, icount, iatom, nraw, ncon + character(len=20) :: start_str + character(len=100) :: lio_dir, line_read + integer, allocatable :: l_of_func(:) + + allocate(l_of_func(MAX_CONTRACT)) + l_of_func = 0 + iostatus = 0 + basis_size = 0 + aux_size = 0 + + if (.not. use_internal) then + ! Read basis from external input basis file. + inquire(file = basis_file, exist = file_exists) + if (.not. file_exists) then + write(*,'(A)') " Error: Basis set file ", trim(basis_file), & + " not found." + iostatus = 1 + return + endif + open(unit = file_uid, file= basis_file, iostat = file_iostat) + read(file_uid, '(A8)') start_str + + do while (start_str .ne. "endbasis") + read(file_uid,*) iatom, nraw, ncon + if (max_f_per_atom .lt. nraw) max_f_per_atom = nraw + if (max_c_per_atom .lt. ncon) max_c_per_atom = ncon + if (ncon .gt. MAX_CONTRACT) then + write(*,'(A)') " Error: Atom has more contractions than the "& + &"maximum allowed (", MAX_CONTRACT,")." + iostatus = 2 + return + endif + + ! Only needs to read angular momenta for each contraction, the rest is + ! skipped. + read(file_uid,*) + read(file_uid,*) (l_of_func(icount), icount=1, ncon) + if (any(atom_Z == iatom)) then + do icount = 1, ncon + basis_size = basis_size & + + ANG_DEG(l_of_func(icount)) * atom_count(iatom) + enddo + endif + do icount = 1, nraw + read(file_uid,*) + enddo + atom_bas_done(iatom) = .true. + + ! Reads auxiliary basis set for an atom doing the same as before. + read(file_uid,*) iatom, nraw, ncon + if (max_f_per_atom .lt. nraw) max_f_per_atom = nraw + if (max_c_per_atom .lt. ncon) max_c_per_atom = ncon + if (ncon .gt. MAX_CONTRACT) then + write(*,'(A)') " Error: Atom has more contractions than the "& + &"maximum allowed (", MAX_CONTRACT,")." + iostatus = 2 + return + endif + + read(file_uid,*) + read(file_uid,*) (l_of_func(icount), icount=1, ncon) + if (any(atom_Z == iatom)) then + do icount = 1, ncon + aux_size = aux_size & + + ANG_DEG(l_of_func(icount)) * atom_count(iatom) + enddo + endif + do icount = 1, nraw + read(file_uid,*) + enddo + atom_fit_done(iatom) = .true. + + read(file_uid, '(A8)') start_str + enddo + close(file_uid) + else + ! Reads from LIO internal basis set files. + call getenv("LIOHOME", lio_dir) + if (lio_dir == "") then + write(*,'(A)') " Error: LIOHOME not set. Cannot use internal basis", & + " files. Either set LIOHOME to your LIO installation", & + " directory or use an external basis set file." + iostatus = 2 + return + endif + basis_file = trim(lio_dir) // "/dat/basis/" // basis_file + fitting_file = trim(lio_dir) // "/dat/basis/fitting/" // fitting_file + + ! Checks file existence. + inquire(file = basis_file, exist = file_exists) + if (.not. file_exists) then + write(*,'(A)') " Error: Basis set file ", trim(basis_file), & + " not found." + iostatus = 1 + return + endif + inquire(file = fitting_file, exist = file_exists) + if (.not. file_exists) then + write(*,'(A)') " Error: Fitting set file ", trim(fitting_file), & + " not found." + iostatus = 1 + return + endif + + ! Reads basis set data. + ! Skips empty lines and those starting with # (comments) + open(unit = file_uid, file= basis_file, iostat = file_iostat) + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (start_str == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + + do while (start_str .ne. "endbasis") + read(file_uid,*) iatom, nraw, ncon + if (max_f_per_atom .lt. nraw) max_f_per_atom = nraw + if (max_c_per_atom .lt. ncon) max_c_per_atom = ncon + if (ncon .gt. MAX_CONTRACT) then + write(*,'(A)') " Error: Atom has more contractions than the "& + &"maximum allowed (", MAX_CONTRACT,")." + iostatus = 2 + return + endif + + ! Only needs to read angular momenta for each contraction, the rest is + ! skipped. + read(file_uid,*) + read(file_uid,*) (l_of_func(icount), icount=1, ncon) + if (any(atom_Z == iatom)) then + do icount = 1, ncon + basis_size = basis_size & + + ANG_DEG(l_of_func(icount)) * atom_count(iatom) + enddo + endif + do icount = 1, nraw + read(file_uid,*) + enddo + atom_bas_done(iatom) = .true. + + ! Skips empty lines or those starting with #. + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (start_str == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + enddo + close(file_uid) + + ! Reads fitting set data. + ! Skips empty lines or those starting with #. + open(unit = file_uid, file= fitting_file, iostat = file_iostat) + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (trim(start_str) == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + + do while (start_str .ne. "endbasis") + read(file_uid,*) iatom, nraw, ncon + if (max_f_per_atom .lt. nraw) max_f_per_atom = nraw + if (max_c_per_atom .lt. ncon) max_c_per_atom = ncon + if (ncon .gt. MAX_CONTRACT) then + write(*,'(A)') " Error: Atom has more contractions than the "& + &"maximum allowed (", MAX_CONTRACT,")." + iostatus = 2 + return + endif + + ! Only needs to read angular momenta for each contraction, the rest is + ! skipped. + read(file_uid,*) + read(file_uid,*) (l_of_func(icount), icount=1, ncon) + if (any(atom_Z == iatom)) then + do icount = 1, ncon + aux_size = aux_size & + + ANG_DEG(l_of_func(icount)) * atom_count(iatom) + enddo + endif + do icount = 1, nraw + read(file_uid,*) + enddo + atom_fit_done(iatom) = .true. + + ! Skips empty lines or those starting with #. + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (trim(start_str) == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + enddo + close(file_uid) + endif + deallocate(l_of_func) +end subroutine basis_set_size + +subroutine check_basis(atom_bas_done, atom_fit_done, n_atoms, atom_Z, iostatus) + implicit none + integer, intent(in) :: n_atoms, atom_Z(n_atoms) + logical, intent(in) :: atom_bas_done(0:120), atom_fit_done(0:120) + integer, intent(out) :: iostatus + + integer :: iatom + character(len=3) :: iatom_name + + iostatus = 0 + do iatom = 1, n_atoms + if ( .not. atom_bas_done(atom_Z(iatom)) ) then + call atom_name(iatom, iatom_name) + write(*,'(A)') " Error: Basis set not found for ", & + trim(iatom_name), "." + iostatus = 1 + return + endif + if ( .not. atom_fit_done(atom_Z(iatom)) ) then + call atom_name(iatom, iatom_name) + write(*,'(A)') " Error: Fitting set not found for ", & + trim(iatom_name), "." + iostatus = 1 + return + endif + enddo +end subroutine check_basis + +subroutine read_basis_external(basis_file, n_funcs, n_fits, n_atoms, normalize,& + max_fun_per_atom, max_con_per_atom, atom_Z, & + coef, expo, coefd, expod, n_cont, n_contd, & + ang_mom_f, ang_mom_fd, atm_of_func, & + atm_of_funcd, min_atm_exp, nns, nnp, nnd, & + nShell, nShelld, iostatus) + use basis_data , only: ANG_DEG + use constants_mod, only: PI32 + implicit none + integer , intent(in) :: max_con_per_atom, max_fun_per_atom, & + n_atoms, atom_Z(n_atoms), n_funcs, n_fits + logical , intent(in) :: normalize + character(len=*), intent(in) :: basis_file + integer , intent(out) :: atm_of_func(n_funcs), atm_of_funcd(n_fits),& + nns(n_atoms), nnp(n_atoms), nnd(n_atoms), & + iostatus, n_cont(n_funcs), n_contd(n_fits),& + ang_mom_f(n_funcs), ang_mom_fd(n_fits), & + nShell(0:3), nShelld(0:3) + double precision, intent(out) :: coef(n_funcs,max_con_per_atom), & + expo(n_funcs,max_con_per_atom), & + coefd(n_fits,max_con_per_atom), & + expod(n_fits,max_con_per_atom), & + min_atm_exp(n_atoms) + + integer :: file_iostat, file_uid = TMP_OPEN_UID + integer :: iatom, nraw, ncon, atom, icont, icount, l2, index, & + n_orig, n_aux + double precision :: min_exp + character(len=20) :: start_str + + logical , allocatable :: basis_done(:), fitting_done(:) + integer , allocatable :: n_cont_func(:), ang_mom(:) + double precision, allocatable :: expo_temp(:), coef_temp(:) + + allocate(n_cont_func(max_con_per_atom) , ang_mom(max_con_per_atom) , & + expo_temp(max_fun_per_atom) , coef_temp(max_fun_per_atom) , & + basis_done(n_atoms) , fitting_done(n_atoms)) + basis_done = .false. + fitting_done = .false. + min_exp = 100000.0D0 + n_orig = 0 + n_aux = 0 + + open(unit = file_uid, file= basis_file, iostat = file_iostat) + read(file_uid,'(A8)') start_str + do while (start_str .ne. 'endbasis') + ! Starts reading the basis set for an atom + ! Reads contraction scheme. The value for p/d/f should not be repeated + ! 3/6/10 times, just set them once. Also reads angular momentum for each + ! of the contractions. + read(file_uid,*) atom, nraw, ncon + read(file_uid,*) n_cont_func(1:ncon) + read(file_uid,*) ang_mom(1:ncon) + do icount = 1, nraw + read(file_uid,*) expo_temp(icount), coef_temp(icount) + if ((expo_temp(icount) .lt. min_exp) .and. any(atom_Z == atom)) & + min_exp = expo_temp(icount) + enddo + + do iatom = 1, n_atoms + if (atom_Z(iatom) .eq. atom .and. (.not. basis_done(iatom))) then + basis_done(iatom) = .true. + min_atm_exp(iatom) = min_exp + + ! These are used for atoms that are near to each other. + nns(iatom) = 0 + nnp(iatom) = 0 + nnd(iatom) = 0 + + do icont = 1, ncon + select case (ang_mom(icont)) + case (0) + nns(iatom) = nns(iatom) + ANG_DEG(ang_mom(icont)) + case (1) + nnp(iatom) = nnp(iatom) + ANG_DEG(ang_mom(icont)) + case (2) + nnd(iatom) = nnd(iatom) + ANG_DEG(ang_mom(icont)) + case default + end select + enddo + + ! Stores the exponents and coefficients, and normalizes them if + ! necessary. + index = 0 + do icont = 1, ncon + nshell(ang_mom(icont)) = nshell(ang_mom(icont)) + & + ANG_DEG(ang_mom(icont)) + do l2 = 1, ANG_DEG(ang_mom(icont)) + n_orig = n_orig +1 + + if (normalize) then + do icount =1, n_cont_func(icont) + index = index +1 + select case (ang_mom(icont)) + case (0) + coef(n_orig, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) / & + PI32) * coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + case (1) + coef(n_orig, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 4.0D0 * expo_temp(index) / & + PI32) * coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + case (2) + coef(n_orig, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 16.0D0 * expo_temp(index)**2/& + PI32) * coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + case default + iostatus = 1 + write(*,'(A)') " ERROR: Basis set contains " , & + "usupported angular momenta. LIO ", & + "uses only s, p and d-type orbitals." + return + end select + enddo + else + do icount =1, n_cont_func(icont) + index = index +1 + coef(n_orig, icount) = coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + enddo + endif + + ! Repeats the index for p, d and f + if (l2 .ne. ANG_DEG(ang_mom(icont))) index = index - & + n_cont_func(icont) + + atm_of_func(n_orig) = iatom + n_cont(n_orig) = n_cont_func(icont) + ang_mom_f(n_orig) = ang_mom(icont) + enddo + enddo + endif + enddo + + ! Starts reading the auxiliary basis for an atom. Same rules as before + ! are applied here. + read(file_uid,*) atom, nraw, ncon + read(file_uid,*) n_cont_func(1:ncon) + read(file_uid,*) ang_mom(1:ncon) + do icount = 1, nraw + read(file_uid,*) expo_temp(icount), coef_temp(icount) + enddo + + do iatom = 1, n_atoms + if (atom_Z(iatom) .eq. atom) then + fitting_done(iatom) = .true. + + index = 0 + do icont = 1, ncon + nshelld(ang_mom(icont)) = nshelld(ang_mom(icont)) + & + ANG_DEG(ang_mom(icont)) + do l2 = 1, ANG_DEG(ang_mom(icont)) + n_aux = n_aux +1 + + if (normalize) then + do icount =1, n_cont_func(icont) + index = index +1 + select case (ang_mom(icont)) + case (0) + coefd(n_aux, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) / & + PI32) * coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + case (1) + coefd(n_aux, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 4.0D0 * expo_temp(index) / & + PI32) * coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + case (2) + coefd(n_aux, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 16.0D0 * expo_temp(index) **2/& + PI32) * coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + case default + iostatus = 1 + write(*,'(A)') " ERROR: Basis set contains " , & + "usupported angular momenta. LIO ", & + "uses only s, p and d-type orbitals." + return + end select + enddo + else + do icount =1, n_cont_func(icont) + index = index +1 + coefd(n_aux, icount) = coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + enddo + endif + + if (l2 .ne. ANG_DEG(ang_mom(icont))) then + index = index - n_cont_func(icont) + endif + + atm_of_funcd(n_aux) = iatom + n_contd(n_aux) = n_cont_func(icont) + ang_mom_fd(n_aux) = ang_mom(icont) + enddo + enddo + endif + enddo + read(file_uid,'(A8)') start_str + enddo +end subroutine read_basis_external + +subroutine read_basis_internal(basis_file, fitting_file, n_funcs, n_fits, & + n_atoms, normalize, max_fun_per_atom, & + max_con_per_atom, atom_Z, coef, expo, coefd, & + expod, n_cont, n_contd, ang_mom_f, ang_mom_fd, & + atm_of_func, atm_of_funcd, min_atm_exp, nns, & + nnp, nnd, nShell, nShelld, iostatus) + use basis_data , only: ANG_DEG + use constants_mod, only: PI32 + implicit none + integer , intent(in) :: max_con_per_atom, max_fun_per_atom, & + n_atoms, atom_Z(n_atoms), n_funcs, n_fits + logical , intent(in) :: normalize + character(len=*), intent(in) :: basis_file, fitting_file + integer , intent(out) :: atm_of_func(n_funcs), atm_of_funcd(n_fits),& + nns(n_atoms), nnp(n_atoms), nnd(n_atoms), & + iostatus, n_cont(n_funcs), n_contd(n_fits),& + ang_mom_f(n_funcs), ang_mom_fd(n_fits), & + nShell(0:3), nShelld(0:3) + double precision, intent(out) :: coef(n_funcs,max_con_per_atom), & + expo(n_funcs,max_con_per_atom), & + coefd(n_fits,max_con_per_atom), & + expod(n_fits,max_con_per_atom), & + min_atm_exp(n_atoms) + + integer :: file_iostat, file_uid = TMP_OPEN_UID + integer :: iatom, nraw, ncon, atom, icont, icount, l2, index, & + n_orig, n_aux + double precision :: min_exp + character(len=20) :: start_str + character(len=100) :: line_read + + logical , allocatable :: basis_done(:), fitting_done(:) + integer , allocatable :: n_cont_func(:), ang_mom(:) + double precision, allocatable :: expo_temp(:), coef_temp(:) + + allocate(n_cont_func(max_con_per_atom +1), ang_mom(max_con_per_atom +1), & + expo_temp(max_fun_per_atom) , coef_temp(max_fun_per_atom) , & + basis_done(n_atoms) , fitting_done(n_atoms)) + basis_done = .false. + fitting_done = .false. + min_exp = 100000.0D0 + n_orig = 0 + n_aux = 0 + + open(unit = file_uid, file= basis_file, iostat = file_iostat) + ! Skips # characters in file. + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (trim(start_str) == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + + do while (start_str .ne. 'endbasis') + ! Starts reading the basis set for an atom + ! Reads contraction scheme. The value for p/d/f should not be repeated + ! 3/6/10 times, just set them once. Also reads angular momentum for each + ! of the contractions. + read(file_uid,*) atom, nraw, ncon + read(file_uid,*) n_cont_func(1:ncon) + read(file_uid,*) ang_mom(1:ncon) + do icount = 1, nraw + read(file_uid,*) expo_temp(icount), coef_temp(icount) + if ((expo_temp(icount) .lt. min_exp) .and. any(atom_Z == atom)) & + min_exp = expo_temp(icount) + enddo + + do iatom = 1, n_atoms + if (atom_Z(iatom) .eq. atom .and. (.not. basis_done(iatom))) then + basis_done(iatom) = .true. + min_atm_exp(iatom) = min_exp + + ! These are used for atoms that are near to each other. + nns(iatom) = 0 + nnp(iatom) = 0 + nnd(iatom) = 0 + + do icont = 1, ncon + select case (ang_mom(icont)) + case (0) + nns(iatom) = nns(iatom) + ANG_DEG(ang_mom(icont)) + case (1) + nnp(iatom) = nnp(iatom) + ANG_DEG(ang_mom(icont)) + case (2) + nnd(iatom) = nnd(iatom) + ANG_DEG(ang_mom(icont)) + case default + end select + enddo + + ! Stores the exponents and coefficients, and normalizes them if + ! necessary. + index = 0 + do icont = 1, ncon + nshell(ang_mom(icont)) = nshell(ang_mom(icont)) + & + ANG_DEG(ang_mom(icont)) + do l2 = 1, ANG_DEG(ang_mom(icont)) + n_orig = n_orig +1 + + if (normalize) then + do icount =1, n_cont_func(icont) + index = index +1 + select case (ang_mom(icont)) + case (0) + coef(n_orig, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) / & + PI32) * coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + case (1) + coef(n_orig, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 4.0D0 * expo_temp(index) / & + PI32) * coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + case (2) + coef(n_orig, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 16.0D0 * expo_temp(index)**2/& + PI32) * coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + case default + iostatus = 1 + write(*,'(A)') " ERROR: Basis set contains " , & + "usupported angular momenta. LIO ", & + "uses only s, p and d-type orbitals." + return + end select + enddo + else + do icount =1, n_cont_func(icont) + index = index +1 + coef(n_orig, icount) = coef_temp(index) + expo(n_orig, icount) = expo_temp(index) + enddo + endif + + ! Repeats the index for p, d and f + if (l2 .ne. ANG_DEG(ang_mom(icont))) index = index - & + n_cont_func(icont) + + atm_of_func(n_orig) = iatom + n_cont(n_orig) = n_cont_func(icont) + ang_mom_f(n_orig) = ang_mom(icont) + enddo + enddo + endif + enddo + + ! Skips # characters in file. + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (trim(start_str) == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + enddo + close(file_uid) + + open(unit = file_uid, file= fitting_file, iostat = file_iostat) + ! Skips # characters in file. + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (trim(start_str) == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + + do while (start_str .ne. 'endbasis') + ! Starts reading the auxiliary basis for an atom. Same rules as before + ! are applied here. + read(file_uid,*) atom, nraw, ncon + read(file_uid,*) n_cont_func(1:ncon) + read(file_uid,*) ang_mom(1:ncon) + do icount = 1, nraw + read(file_uid,*) expo_temp(icount), coef_temp(icount) + enddo + + do iatom = 1, n_atoms + if (atom_Z(iatom) .eq. atom) then + fitting_done(iatom) = .true. + + index = 0 + do icont = 1, ncon + nshelld(ang_mom(icont)) = nshelld(ang_mom(icont)) + & + ANG_DEG(ang_mom(icont)) + do l2 = 1, ANG_DEG(ang_mom(icont)) + n_aux = n_aux +1 + + if (normalize) then + do icount =1, n_cont_func(icont) + index = index +1 + select case (ang_mom(icont)) + case (0) + coefd(n_aux, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) / & + PI32) * coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + case (1) + coefd(n_aux, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 4.0D0 * expo_temp(index) / & + PI32) * coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + case (2) + coefd(n_aux, icount) = dsqrt( dsqrt(8.0D0 * & + (expo_temp(index)) ** 3 ) * & + 16.0D0 * expo_temp(index)**2/& + PI32) * coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + case default + iostatus = 1 + write(*,'(A)') " ERROR: Basis set contains " , & + "usupported angular momenta. LIO ", & + "uses only s, p and d-type orbitals." + return + end select + enddo + else + do icount =1, n_cont_func(icont) + index = index +1 + coefd(n_aux, icount) = coef_temp(index) + expod(n_aux, icount) = expo_temp(index) + enddo + endif + + if (l2 .ne. ANG_DEG(ang_mom(icont))) then + index = index - n_cont_func(icont) + endif + + atm_of_funcd(n_aux) = iatom + n_contd(n_aux) = n_cont_func(icont) + ang_mom_fd(n_aux) = ang_mom(icont) + enddo + enddo + endif + enddo + + ! Skips # characters in file. + line_read = "" + do while (line_read == "") + read(file_uid,*) line_read + read(line_read, '(A1)') start_str + if (trim(start_str) == "#") line_read = "" + enddo + read(line_read, '(A8)') start_str + enddo + +end subroutine read_basis_internal + +subroutine reorder_basis(expon, coeff, atom_of_funct, n_cont, mixed_index, & + basis_size, max_cont, l_of_funct, n_shell) + implicit none + integer , intent(in) :: basis_size, max_cont, n_shell(0:3), & + l_of_funct(basis_size) + integer , intent(inout) :: atom_of_funct(basis_size), & + n_cont(basis_size) , & + mixed_index(basis_size) + double precision, intent(inout) :: expon(basis_size, max_cont), & + coeff(basis_size, max_cont) + + double precision, allocatable :: expo_t(:,:), coef_t(:,:) + integer , allocatable :: atom_of_funct_t(:), n_cont_t(:) + integer :: ifunct, s_index, p_index, d_index + + allocate(expo_t(basis_size, max_cont), coef_t(basis_size, max_cont), & + atom_of_funct_t(basis_size) , n_cont_t(basis_size)) + + s_index = 1 + p_index = 1 + n_shell(0) + d_index = 1 + n_shell(0) + n_shell(1) + + do ifunct = 1, basis_size + select case (l_of_funct(ifunct)) + case (0) ! s function + atom_of_funct_t(s_index) = atom_of_funct(ifunct) + mixed_index(s_index) = ifunct + n_cont_t(s_index) = n_cont(ifunct) + expo_t(s_index,:) = expon(ifunct,:) + coef_t(s_index,:) = coeff(ifunct,:) + + s_index = s_index +1 + case (1) ! p functions + atom_of_funct_t(p_index) = atom_of_funct(ifunct) + mixed_index(p_index) = ifunct + n_cont_t(p_index) = n_cont(ifunct) + expo_t(p_index,:) = expon(ifunct,:) + coef_t(p_index,:) = coeff(ifunct,:) + + p_index = p_index +1 + case (2) ! d functions + atom_of_funct_t(d_index) = atom_of_funct(ifunct) + mixed_index(d_index) = ifunct + n_cont_t(d_index) = n_cont(ifunct) + expo_t(d_index,:) = expon(ifunct,:) + coef_t(d_index,:) = coeff(ifunct,:) + + d_index = d_index +1 + case default + end select + enddo + + n_cont = n_cont_t + atom_of_funct = atom_of_funct_t + expon = expo_t + coeff = coef_t + deallocate(expo_t, coef_t, atom_of_funct_t, n_cont_t) +end subroutine reorder_basis +!############################################################################### +! TAKE THIS TO OTHER MODULE +subroutine atom_name(atom_Z, symb) + ! Takes atomic number Z and translates it to its name. + implicit none + integer , intent(in) :: atom_Z + character(LEN=3), intent(out) :: symb + + character(LEN=3) :: name(118) + name = (/'H ','HE ','LI ','BE ','B ','C ','N ','O ','F ','NE ','NA ', & + 'MG ','AL ','SI ','P ','S ','CL ','AR ','K ','CA ','SC ','TI ', & + 'V ','CR ','MN ','FE ','CO ','NI ','CU ','ZN ','GA ','GE ','AS ', & + 'SE ','BR ','KR ','RB ','SR ','Y ','ZR ','NB ','MO ','TC ','RU ', & + 'RH ','PD ','AG ','CD ','IN ','SN ','SB ','TE ','I ','XE ','CS ', & + 'BA ','LA ','CE ','PR ','ND ','PM ','SM ','EU ','GD ','TB ','DY ', & + 'HO ','ER ','TM ','YB ','LU ','HF ','TA ','W ','RE ','OS ','IR ', & + 'PT ','AU ','HG ','TL ','PB ','BI ','PO ','AT ','RN ','FR ','RA ', & + 'AC ','TH ','PA ','U ','NP ','PU ','AM ','CM ','BK ','CF ','ES ', & + 'FM ','MD ','NO ','LR ','RF ','DB ','SG ','BH ','HS ','MT ','DS ', & + 'UUU','UUB','UUT','UUQ','UUP','UUH','UUS','UUO'/) + symb = name(atom_Z) + +end subroutine atom_name + +end module basis_subs !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! diff --git a/lioamber/basis_subs.f90 b/lioamber/basis_subs.f90 deleted file mode 100644 index 3d7f1969b..000000000 --- a/lioamber/basis_subs.f90 +++ /dev/null @@ -1,156 +0,0 @@ -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! -module basis_subs - implicit none - contains -!------------------------------------------------------------------------------! -! -! -! -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! -subroutine basis_data_set(ns,np,nd,orba,orbc,ge,gc) -!------------------------------------------------------------------------------! - use basis_data, only: basis_size, basis_size_s, basis_size_p, basis_size_d & - &, maximum_contractions, orbital_contractions & - &, parent_atom, angular_momentum, gauss_expo, gauss_coef - implicit none - integer,intent(in) :: ns,np,nd - integer,intent(in), dimension(:) :: orba - integer,intent(in), dimension(:) :: orbc - real*8,intent(in) , dimension(:,:) :: ge - real*8,intent(in) , dimension(:,:) :: gc - integer :: nn - - - basis_size_s=ns - basis_size_p=np - basis_size_d=nd -! basis_size_f=nf - basis_size=ns+np+nd - maximum_contractions=max(size(ge,2),size(gc,2)) - - if (allocated(parent_atom)) deallocate(parent_atom) - if (allocated(orbital_contractions)) deallocate(orbital_contractions) - if (allocated(angular_momentum)) deallocate(angular_momentum) - if (allocated(gauss_expo)) deallocate(gauss_expo) - if (allocated(gauss_coef)) deallocate(gauss_coef) - - allocate( parent_atom(basis_size) ) - allocate( orbital_contractions(basis_size) ) - allocate( angular_momentum(3, basis_size) ) - allocate( gauss_expo(maximum_contractions, basis_size) ) - allocate( gauss_coef(maximum_contractions, basis_size) ) - - - do nn=1,ns - parent_atom(nn+0) = orba(nn) - orbital_contractions(nn+0) = orbc(nn) - - angular_momentum(:,nn) = 0 - gauss_expo(:,nn+0) = ge(nn,:) - gauss_coef(:,nn+0) = gc(nn,:) - enddo - - do nn=ns+1,ns+np,3 - parent_atom(nn+2) = orba(nn) - parent_atom(nn+1) = orba(nn) - parent_atom(nn+0) = orba(nn) - orbital_contractions(nn+2) = orbc(nn) - orbital_contractions(nn+1) = orbc(nn) - orbital_contractions(nn+0) = orbc(nn) - - angular_momentum(:,nn+2) = 0 - angular_momentum(:,nn+1) = 0 - angular_momentum(:,nn+0) = 0 - angular_momentum(1,nn+0) = 1 ! px - angular_momentum(2,nn+1) = 1 ! py - angular_momentum(3,nn+2) = 1 ! pz - - gauss_expo(:,nn+2) = ge(nn,:) - gauss_expo(:,nn+1) = ge(nn,:) - gauss_expo(:,nn+0) = ge(nn,:) - - gauss_coef(:,nn+2) = gc(nn,:) - gauss_coef(:,nn+1) = gc(nn,:) - gauss_coef(:,nn+0) = gc(nn,:) - enddo - - do nn=ns+np+1,ns+np+nd,6 - parent_atom(nn+5) = orba(nn) - parent_atom(nn+4) = orba(nn) - parent_atom(nn+3) = orba(nn) - parent_atom(nn+2) = orba(nn) - parent_atom(nn+1) = orba(nn) - parent_atom(nn+0) = orba(nn) - orbital_contractions(nn+5) = orbc(nn) - orbital_contractions(nn+4) = orbc(nn) - orbital_contractions(nn+3) = orbc(nn) - orbital_contractions(nn+2) = orbc(nn) - orbital_contractions(nn+1) = orbc(nn) - orbital_contractions(nn+0) = orbc(nn) - - angular_momentum(:,nn+5) = 0 - angular_momentum(:,nn+4) = 0 - angular_momentum(:,nn+3) = 0 - angular_momentum(:,nn+2) = 0 - angular_momentum(:,nn+1) = 0 - angular_momentum(:,nn+0) = 0 - angular_momentum(1,nn+0) = 2 ! dxx (x) - angular_momentum(1,nn+1) = 1 ! dxy (x) - angular_momentum(2,nn+1) = 1 ! dxy (y) - angular_momentum(2,nn+2) = 2 ! dyy (y) - angular_momentum(1,nn+3) = 1 ! dxz (x) - angular_momentum(3,nn+3) = 1 ! dxz (z) - angular_momentum(2,nn+4) = 1 ! dyz (y) - angular_momentum(3,nn+4) = 1 ! dyz (z) - angular_momentum(3,nn+5) = 2 ! dzz (z) - - gauss_expo(:,nn+5) = ge(nn,:) - gauss_expo(:,nn+4) = ge(nn,:) - gauss_expo(:,nn+3) = ge(nn,:) - gauss_expo(:,nn+2) = ge(nn,:) - gauss_expo(:,nn+1) = ge(nn,:) - gauss_expo(:,nn+0) = ge(nn,:) - - gauss_coef(:,nn+5) = gc(nn,:)/SQRT(3.0) - gauss_coef(:,nn+4) = gc(nn,:) - gauss_coef(:,nn+3) = gc(nn,:) - gauss_coef(:,nn+2) = gc(nn,:)/SQRT(3.0) - gauss_coef(:,nn+1) = gc(nn,:) - gauss_coef(:,nn+0) = gc(nn,:)/SQRT(3.0) - enddo - -end subroutine -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! -subroutine basis_data_norm( Isize, Icont, gcoefs ) -!------------------------------------------------------------------------------! - use maskrmm, only: rmmget_fock - use faint_cpu, only: int1 - use basis_data, only: gauss_coef - use garcha_mod, only: RMM,Nuc,a,c,d,r,Iz,NORM,natom,M,Md,ncont,nshell,ntatom - - implicit none - integer, intent(in) :: Isize - integer, intent(in) :: Icont - double precision , intent(inout) :: gcoefs(Isize, Icont) - - integer :: kk, MM, M5, M11 - double precision :: scratch_energy - double precision, allocatable :: Smat(:,:) - -! call g2g_timer_start('RMMcalc1') - allocate( Smat(isize, isize) ) - MM = M*(M+1)/2 ; M5 = 1 + MM*2; M11 = M5+MM+Md*(Md+1) - call int1(scratch_energy, RMM(M5:M5+MM), RMM(M11:M11+MM), Smat, d, r, Iz, & - natom, ntatom) - do kk = 1, isize - gcoefs(kk,:) = gcoefs(kk,:) / sqrt( Smat(kk,kk) ) - gauss_coef(:,kk) = gauss_coef(:,kk) / sqrt( Smat(kk,kk) ) - enddo - deallocate( Smat ) -! call g2g_timer_stop('RMMcalc1') - -end subroutine basis_data_norm - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! -end module -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! diff --git a/lioamber/drive.f90 b/lioamber/drive.f90 index e8c3995c0..780cfe648 100644 --- a/lioamber/drive.f90 +++ b/lioamber/drive.f90 @@ -8,7 +8,7 @@ ! !------------------------------------------------------------------------------! - SUBROUTINE drive(ng2,ngDyn,ngdDyn) +subroutine drive(ng2,ngDyn,ngdDyn) USE garcha_mod, ONLY : a,c, basis, done, done_fit, natomc, nnps, & nnpp, nnpd, nns, nnp, nnd, atmin, jatc, ncf, lt, at, ct, nnat, nshell, & nuc, ncont, nlb, nshelld, cd, ad, Nucd, ncontd, nld, Nucx, indexii, & @@ -69,6 +69,7 @@ SUBROUTINE drive(ng2,ngDyn,ngdDyn) ! ! NORM true , expansion in normalized gaussians, so normalization factor ! included in coefficients + nopt=0 ! ! calls generator of table for incomplete gamma functions @@ -80,820 +81,9 @@ SUBROUTINE drive(ng2,ngDyn,ngdDyn) !----------------------------------------------------------------------- ! reads input file - - if (.not.int_basis) then - inquire(file=basis,exist=exists) - if (.not.exists) then - write(*,*) 'ERROR CANNOT FIND INPUT FILE ON UNIT 1',basis - stop - else -! name of output file - ikk=1 - do while (basis(ikk:ikk).ne.' ') - ikk=ikk+1 - enddo -! - ikk=ikk-1 -! name2=basis(1:ikk)//'.out' -! - open(unit=1,file=basis,iostat=ios) -! open(unit=2,file=output) - endif - endif - - if (writexyz) open(unit=18,file=fcoord) - if ((mulliken).or.(td_do_pop.gt.0)) open(unit=85,file=fmulliken) - if (restart_freq.gt.0) open(unit=88,file=frestart) - -!------------------------------------------------------- - do i=1,natom - done(i)=.false. - done_fit(i)=.false. - enddo - -!c ------------------------------------------------------------- -!c for each kind of atom, given by Z -!c -!c Basis set for MO expansion -!c -!c reads whatis: Gaussian for the time being, -!c later stdbas will be included -!c reads iatom, atomic number , nraw number of primitive gaussians -!c , counting 3 p 6 d .., and ncon number of contractions counted in the -!c same way -!c if coefficients correspond to normalized gaussians -!c normalization factor will multiply coefficients, so in the program -!c expressions for non-normalized gaussians are used. -!c an option is included for using non-normalized gaussians -!c !The contractions don't need to be normalized, it will be done -!c automatically -!c -!c - if (.not.int_basis) then - read(1,100) whatis - endif - - No=0 - Nd=0 - NBAS=0 - M=0 - Md=0 - - allocate(natomc(natom),nnps(natom),nnpp(natom),nnp(natom)) - allocate(nnpd(natom),nns(natom),nnd(natom),atmin(natom)) - allocate(jatc(natom,natom)) - - do i=1,natom - natomc(i)=0 - enddo -!c------------------------------------------------------------------------- -!c BASIS SETS ------------------------------------------------------------ -!c------------------------------------------------------------------------- - if (.not.int_basis) then - do 25 while (whatis.ne.'endbasis') -!c - NBAS=NBAS+1 -!c signals if a basis set was not used - used=.false. - atmint=100000. - read(1,*) iatom,nraw,ncon -!c write(2,600) iatom,nraw,ncon -!c -!c reads contraction scheme. The value for p,d ,f should not be repeated -!c 3 ,6 , 10 ..... times -!c reads also angular momentum for each of the contractions -!c 0 for s 1 for p, etc -!c - read(1,*) (ncf(i),i=1,ncon) -!c write(2,*) (ncf(i),i=1,ncon) - read(1,*) (lt(i),i=1,ncon) -!c write(2,*) (lt(i),i=1,ncon) -!c -!c loop over all primitives, no repeating p, d - do i=1,nraw - read(1,*) at(i),ct(i) - if(at(i).lt.atmint) atmint=at(i) - enddo -!c - do j=1,natom - if(Iz(j).eq.iatom.and.(.not.done(j))) then - nnat(NBAS)=nnat(NBAS)+1 - done(j)=.true. - used=.true. - atmin(j)=atmint - -!c cosas que puso nano para "atomos cerca" - - nns(j)=0 - nnp(j)=0 - nnd(j)=0 - - do kkk=1,ncon - if (lt(kkk).eq.0) nns(j)=nns(j)+Num(lt(kkk)) - if (lt(kkk).eq.1) nnp(j)=nnp(j)+Num(lt(kkk)) - if (lt(kkk).eq.2) nnd(j)=nnd(j)+Num(lt(kkk)) - enddo - -!c write(*,*) 'nns y etc',nns(j),nnp(j),nnd(j),nnps(j) -!c > ,nnpp(j),nnpd(j) - -!c =====>>>>>> M stores # of contractions <<<<<<=========== - index=0 - do k=1,ncon -!c - M=M+Num(lt(k)) - -!c nshell gives the # of functions s, p, d etc - nshell(lt(k))=nshell(lt(k))+Num(lt(k)) -!c - do l2=1,Num(lt(k)) - No=No+1 -!c -!c normalization -!c - if (NORM) then - do l=1,ncf(k) -!c - index=index+1 - - if(lt(k).eq.0) then -!c - xnorm=sqrt((2.D0*at(index)/pi)**3) - xnorm=sqrt(xnorm) - c(No,l)=ct(index)*xnorm - a(No,l)=at(index) - elseif(lt(k).eq.1) then -!c - xnorm=sqrt((2.D0*at(index)/pi)**3)*4.D0*at(index) - xnorm=sqrt(xnorm) - c(No,l)=ct(index)*xnorm - a(No,l)=at(index) - elseif (lt(k).eq.2) then -!c - xnorm=sqrt((2.D0*at(index)/pi)**3)*(4.D0*at(index))**2 - xnorm=sqrt(xnorm) - c(No,l)=ct(index)*xnorm - a(No,l)=at(index) - endif - enddo - else -!c no normalization case - do l=1,ncf(k) - index=index+1 - c(No,l)=ct(index) - a(No,l)=at(index) - enddo - endif - -!c repeat the index only for p,d and f - if (l2.ne.Num(lt(k))) then - index=index-ncf(k) - endif -!c - Nuc(No)=j - ncont(No)=ncf(k) - nlb(No)=lt(k) - enddo - enddo - endif - enddo -!c - if ( (.not.used) .and. (verbose.gt.4) ) then - write(*,200) iatom - endif -!c -!c Exactly the same should be repeated for charge density -!c and exchange correlation potential basis sets -!c -!c CHARGE DENSITY -------------------------------------------------- -!c - read(1,*) iatom,nraw,ncon -!c -!c reads contraction scheme. The value for p,d ,f should not be repeated -!c 3 ,6 , 10 ..... times. Reads also angular type , -!c 0 for s , 1 for p etc - read(1,*) (ncf(i),i=1,ncon) - read(1,*) (lt(i),i=1,ncon) -!c -!c loop over all primitives, repeating p, d - do i=1,nraw - read(1,*) at(i),ct(i) - enddo -!c - do 45 j=1,natom - if (Iz(j).eq.iatom) then - done_fit(j)=.true. -!c -!c Mdd stores # of contractions in final basis, counting all possibilities -!c for p , d etc -!c - index=0 - do 46 k=1,ncon -!c - Md=Md+Num(lt(k)) - nshelld(lt(k))=nshelld(lt(k))+Num(lt(k)) -!c - do 46 l2=1,Num(lt(k)) -!c - Nd=Nd+1 -!c - if (NORM) then - do 47 l=1,ncf(k) - index=index+1 -!c - goto (71,81,91) lt(k)+1 -!c - 71 xnorm=sqrt((2.D0*at(index)/pi)**3) - xnorm=sqrt(xnorm) - cd(Nd,l)=ct(index)*xnorm - ad(Nd,l)=at(index) -!c ad(Nd,l)=2.D0*at(index) - goto 47 -!c - 81 xnorm=sqrt((2.D0*at(index)/pi)**3)*4.D0*at(index) - xnorm=sqrt(xnorm) - cd(Nd,l)=ct(index)*xnorm -!c ad(Nd,l)=2.D0*at(index) - ad(Nd,l)=at(index) - goto 47 -!c - 91 xnorm=sqrt((2.D0*at(index)/pi)**3)*(4.D0*at(index))**2 - xnorm=sqrt(xnorm) - cd(Nd,l)=ct(index)*xnorm -!c ad(Nd,l)=2.D0*at(index) - ad(Nd,l)=at(index) - goto 47 -!c - 47 continue - else -!c -!c no normalization case -!c - do l=1,ncf(k) - index=index+1 -!c - cd(Nd,l)=ct(index) -!c ad(Nd,l)=2.D0*at(index) - ad(Nd,l)=at(index) - enddo - endif -!c -!c repeat the index only for p,d and f, criterium l2>>>>> M stores # of contractions <<<<<<=========== - index=0 - do k=1,ncon -!c - M=M+Num(lt(k)) - -!c nshell gives the # of functions s, p, d etc - nshell(lt(k))=nshell(lt(k))+Num(lt(k)) -!c - do l2=1,Num(lt(k)) - No=No+1 -!c -!c normalization -!c - if (NORM) then - do l=1,ncf(k) -!c - index=index+1 - - if(lt(k).eq.0) then -!c - xnorm=sqrt((2.D0*at(index)/pi)**3) - xnorm=sqrt(xnorm) - c(No,l)=ct(index)*xnorm - a(No,l)=at(index) - elseif(lt(k).eq.1) then -!c - xnorm=sqrt((2.D0*at(index)/pi)**3)*4.D0*at(index) - xnorm=sqrt(xnorm) - c(No,l)=ct(index)*xnorm - a(No,l)=at(index) - elseif (lt(k).eq.2) then -!c - xnorm=sqrt((2.D0*at(index)/pi)**3)*(4.D0*at(index))**2 - xnorm=sqrt(xnorm) - c(No,l)=ct(index)*xnorm - a(No,l)=at(index) - else - write(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - write(*,*) "The basis set ",trim(basis_set), & - " contains f (or higher) functions, and lio does not currently", & - " support them. Choose another basis set." - write(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - stop - endif - enddo - else -!c no normalization case - do l=1,ncf(k) - index=index+1 - c(No,l)=ct(index) - a(No,l)=at(index) - enddo - endif - -!c repeat the index only for p,d and f - if (l2.ne.Num(lt(k))) then - index=index-ncf(k) - endif -!c - Nuc(No)=j - ncont(No)=ncf(k) - nlb(No)=lt(k) - enddo - enddo - endif - enddo -!c - if (.not.used.and.(verbose.gt.4).and. .not.omit_bas) then - write(*,200) iatom - endif - - ! skip over all commented and blank lines - inp_line = "" - do while (inp_line == "") - read(1,*) inp_line - read(inp_line,'(a1)') inp_char - if (inp_char == "#" .or. inp_line == "") then - inp_line = "" - endif - enddo - read(inp_line,100) whatis - -26 enddo - - close(1) - - open(unit=1,file=fit_basis_file,iostat=ios) - - ! skip over all commented and blank lines - inp_line = "" - do while (inp_line == "") - read(1,*) inp_line - read(inp_line,'(a1)') inp_char - if (inp_char == "#" .or. inp_line == "") then - inp_line = "" - endif - enddo - - read(inp_line,100) whatis -!c -!c Exactly the same should be repeated for charge density -!c and exchange correlation potential basis sets -!c -!c CHARGE DENSITY -------------------------------------------------- -!c - ! NOTE: we don't check right now for total primitives in the cd basis - ! being less than nng...probably not a problem... - do 27 while (whatis.ne.'endbasis') - read(1,*) iatom,nraw,ncon -!c write(2,*) iatom,nraw,ncon -!c -!c reads contraction scheme. The value for p,d ,f should not be repeated -!c 3 ,6 , 10 ..... times. Reads also angular type , -!c 0 for s , 1 for p etc - read(1,*) (ncf(i),i=1,ncon) - read(1,*) (lt(i),i=1,ncon) -!c -!c loop over all primitives, repeating p, d - do i=1,nraw - read(1,*) at(i),ct(i) - enddo -!c - do 48 j=1,natom - if (Iz(j).eq.iatom) then - done_fit(j)=.true. -!c -!c Mdd stores # of contractions in final basis, counting all possibilities -!c for p , d etc -!c - index=0 - do 49 k=1,ncon -!c - Md=Md+Num(lt(k)) - nshelld(lt(k))=nshelld(lt(k))+Num(lt(k)) -!c - do 49 l2=1,Num(lt(k)) -!c - Nd=Nd+1 -!c - if (NORM) then - do 50 l=1,ncf(k) - index=index+1 -!c - goto (72,82,92) lt(k)+1 - write(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - write(*,*) "The basis set ",trim(fitting_set), & - " contains f (or higher) functions, and lio does not currently", & - " support them. Choose another basis set" - write(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - stop -!c - 72 xnorm=sqrt((2.D0*at(index)/pi)**3) - xnorm=sqrt(xnorm) - cd(Nd,l)=ct(index)*xnorm - ad(Nd,l)=at(index) - goto 50 -!c - 82 xnorm=sqrt((2.D0*at(index)/pi)**3)*4.D0*at(index) - xnorm=sqrt(xnorm) - cd(Nd,l)=ct(index)*xnorm -!c ad(Nd,l)=2.D0*at(index) - ad(Nd,l)=at(index) - goto 50 -!c - 92 xnorm=sqrt((2.D0*at(index)/pi)**3)*(4.D0*at(index))**2 - xnorm=sqrt(xnorm) - cd(Nd,l)=ct(index)*xnorm -!c ad(Nd,l)=2.D0*at(index) - ad(Nd,l)=at(index) - goto 50 -!c - 50 continue - else -!c -!c no normalization case -!c - do l=1,ncf(k) - index=index+1 -!c - cd(Nd,l)=ct(index) -!c ad(Nd,l)=2.D0*at(index) - ad(Nd,l)=at(index) - enddo - endif -!c -!c repeat the index only for p,d and f, criterium l2>>>>> M stores # of contractions <<<<<<=========== - END DO - END IF - END DO - - inp_line = "" - DO WHILE (inp_line == "") ! skip over all commented and blank lines - READ(1,*) inp_line - READ(inp_line,'(a1)') inp_char - IF (inp_char == "#" .OR. inp_line == "") THEN - inp_line = "" - END IF - END DO - READ(inp_line,100) whatis -26 END DO - CLOSE(1) - -!c -!c Exactly the same should be repeated for charge density -!c -!c CHARGE DENSITY -------------------------------------------------- - - OPEN(UNIT=1,FILE=fit_basis_file,IOSTAT=ios) - - inp_line = "" - DO WHILE (inp_line == "") ! skip over all commented and blank lines - READ(1,*) inp_line - READ(inp_line,'(a1)') inp_char - IF (inp_char == "#" .OR. inp_line == "") THEN - inp_line = "" - END IF - END DO -! - READ(inp_line,100) whatis - - DO 27 WHILE (whatis .NE. 'endbasis') - READ(1,*) iatom,nraw,ncon - - IF (max_func .LT. nraw) max_func=nraw - IF (ncon .GE. nng) call reallocate_ncf_lt(ncon) !agregada para variables dinamicas - - READ(1,*) (ncf(i),i=1,ncon) - READ(1,*) (lt(i),i=1,ncon) - DO i=1,nraw - READ(1,*) !at(i),ct(i) - END DO -!c - DO 48 j=1,natom - IF (Iz(j).EQ.iatom) THEN - done_fit(j)=.true. - DO k=1,ncon - Md=Md+Num(lt(k)) - END DO - END IF - 48 END DO - - inp_line = "" - DO WHILE (inp_line == "") ! skip over all commented and blank lines - READ(1,*) inp_line - READ(inp_line,'(a1)') inp_char - IF (inp_char == "#" .or. inp_line == "") THEN - inp_line = "" - END IF - END DO -! - READ(inp_line,100) whatis - 27 END DO - CLOSE(1) - END IF - - -!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - basis_check=.true. - DO i=1,natom ! chekeo de lectura de base para todos los atomos, Nick - IF (.NOT. done(i)) THEN - CALL asignacion(Iz(i),simb) - WRITE(*,*) simb," havent got a basis set, check basis file" - basis_check=.false. - END IF - END DO -!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - DO i=1,natom ! chekeo de lectura de base auxiliar para todos los atomos, Nick - IF (.NOT. done_fit(i)) THEN - CALL asignacion(Iz(i),simb) - WRITE(*,*) simb," havent got a fitting basis set, check", & - " auxiliar basis file" - basis_check=.false. - END IF - END DO - IF (.NOT. basis_check) STOP -!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ngnu=M - ngdnu=Md - ALLOCATE (at(max_func), ct(max_func)) - - 100 FORMAT (A8) - RETURN - END SUBROUTINE DIMdrive - - - SUBROUTINE reallocate_ncf_lt(Dime) - USE garcha_mod, ONLY: ncf, lt, nng - IMPLICIT NONE - INTEGER, INTENT(IN) :: Dime - nng=Dime+1 - DEALLOCATE (ncf, lt) - ALLOCATE (ncf(nng), lt(nng)) - ENDSUBROUTINE reallocate_ncf_lt +end subroutine drive subroutine get_nco(atom_Z, n_atoms, n_orbitals, n_unpaired, charge, open_shell) use ghost_atoms_subs, only: adjust_ghost_charge diff --git a/lioamber/ehrensubs/calc_forceDS_dds.f90 b/lioamber/ehrensubs/calc_forceDS_dds.f90 index a14c1997f..b2c9e8a3a 100644 --- a/lioamber/ehrensubs/calc_forceDS_dds.f90 +++ b/lioamber/ehrensubs/calc_forceDS_dds.f90 @@ -5,7 +5,9 @@ subroutine calc_forceDS_dds( natoms, nbasis, pos, vel, Mat0, fterm ) ! DESCRIPTION ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! - use basis_data + use basis_data, only: angular_momentum=>ang_mom_ehren, gauss_coef=>c_ehren, & + gauss_expo=>a_ehren + use garcha_mod, only: parent_atom=>nuc, orbital_contractions=>nCont implicit none integer,intent(in) :: natoms ! Number of atoms integer,intent(in) :: nbasis ! Number of basis diff --git a/lioamber/ehrensubs/calc_forceDS_dss.f90 b/lioamber/ehrensubs/calc_forceDS_dss.f90 index fdc6feb3c..a2a0f5805 100644 --- a/lioamber/ehrensubs/calc_forceDS_dss.f90 +++ b/lioamber/ehrensubs/calc_forceDS_dss.f90 @@ -5,7 +5,9 @@ subroutine calc_forceDS_dss(natoms,nbasis,pos,vel,Mat0,MatB,fterm) ! DESCRIPTION ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! - use basis_data + use basis_data, only: angular_momentum=>ang_mom_ehren, gauss_coef=>c_ehren, & + gauss_expo=>a_ehren + use garcha_mod, only: parent_atom=>nuc, orbital_contractions=>nCont implicit none integer,intent(in) :: natoms ! Number of atoms integer,intent(in) :: nbasis ! Number of basis diff --git a/lioamber/func.f b/lioamber/func.f index 60cad8e70..e677a6185 100644 --- a/lioamber/func.f +++ b/lioamber/func.f @@ -10,50 +10,50 @@ c program c 11 March 1992 ----------------------------------------------- c - FUNCTION FUNCT(N,T) + FUNCTION FUNCT(N,T) USE garcha_mod - IMPLICIT REAL*8 (A-H,O-Z) + IMPLICIT REAL*8 (A-H,O-Z) ! COMMON /TABLE/ STR(880,0:21),FAC(0:16) C if (T.lt.0.0D0) then write(*,*) 'Problems',T T=abs(T) endif -C - IF (T.LE.43.975D0) THEN - IT = 20.0D0 * (T + 0.025D0) - TI = DFLOAT(IT) - IT=IT + 1 - DELT = T - 0.05D0 * TI - DELT3 = DELT * 0.333333333333333D0 - DELT4 = 0.25D0 * DELT - DELT2 = DELT4 + DELT4 - DELT5 = 0.20D0 * DELT -c -c -c - TF0 = STR(IT,N) - TF1 = STR(IT,N+1) - TF2 = STR(IT,N+2) - TF3 = STR(IT,N+3) - TF4 = STR(IT,N+4) - TF5 = STR(IT,N+5) -c - FCAP=TF0-DELT*( TF1-DELT2*(TF2-DELT3*(TF3-DELT4* - > (TF4-DELT5*TF5)))) - FUNCT = FCAP - RETURN -C - +C + IF (T.LE.43.975D0) THEN + IT = 20.0D0 * (T + 0.025D0) + TI = DFLOAT(IT) + IT=IT + 1 + DELT = T - 0.05D0 * TI + DELT3 = DELT * 0.333333333333333D0 + DELT4 = 0.25D0 * DELT + DELT2 = DELT4 + DELT4 + DELT5 = 0.20D0 * DELT +c +c +c + TF0 = STR(IT,N) + TF1 = STR(IT,N+1) + TF2 = STR(IT,N+2) + TF3 = STR(IT,N+3) + TF4 = STR(IT,N+4) + TF5 = STR(IT,N+5) +c + FCAP=TF0-DELT*( TF1-DELT2*(TF2-DELT3*(TF3-DELT4* + > (TF4-DELT5*TF5)))) + FUNCT = FCAP + RETURN +C + ELSE FUNCT=FAC(N)*1.D0/(T**N*dsqrt(T)) - + ENDIF c - END - + END -C ---------------------------------------------- + +C ---------------------------------------------- c subroutine for generating tables, used later on in the Taylor expansion c for the incomplete Gamma functions c----------------------------------------------- @@ -61,6 +61,7 @@ SUBROUTINE GENERF USE garcha_mod IMPLICIT REAL*8 (A-H,O-Z) parameter (sqpi=1.77245385090551588D0) + integer icount ! COMMON /TABLE/ STR(880,0:21),FAC(0:16) c c loop over T values in the table ( 0. to 43.95 , interval 0.05) @@ -74,10 +75,10 @@ SUBROUTINE GENERF c loop over order of incomple Gamma functions ( 0 to 21, the ones c necessary for evaluating orders 0-16) c - DO 99 M=21,0,-1 - W=2.D0*DFLOAT(M)+1.D0 - STR(I,M) = (Y + U * F)/W - F=STR(I,M) + DO 99 icount=21,0,-1 + W=2.D0*DFLOAT(icount)+1.D0 + STR(I,icount) = (Y + U * F)/W + F=STR(I,icount) 99 CONTINUE c 799 CONTINUE @@ -94,6 +95,7 @@ SUBROUTINE GENERF SUBROUTINE GENERFS IMPLICIT REAL*4 (A-H,O-Z) parameter (sqpi=1.772453850905E0) + integer icount COMMON /TABLES/ STRS(880,0:21),FACS(0:16) c c loop over T values in the table ( 0. to 43.95 , interval 0.05) @@ -107,10 +109,10 @@ SUBROUTINE GENERFS c loop over order of incomple Gamma functions ( 0 to 21, the ones c necessary for evaluating orders 0-16) c - DO 99 M=21,0,-1 - W=2.E0*FLOAT(M)+1.E0 - STRS(I,M) = (Y + U * F)/W - F=STRS(I,M) + DO 99 icount=21,0,-1 + W=2.E0*FLOAT(icount)+1.E0 + STRS(I,icount) = (Y + U * F)/W + F=STRS(I,icount) 99 CONTINUE c 799 CONTINUE @@ -140,10 +142,10 @@ FUNCTION FMCH(M,X) TERM = TERM * X / (A + DFLOAT(I-1)) PTLSUM = PTLSUM + TERM IF(TERM/PTLSUM.LT.1.0D-12)GO TO 12 - 11 CONTINUE - STOP + 11 CONTINUE + STOP 12 FMCH = 0.5D0 * PTLSUM * Y - RETURN + RETURN c 20 A = DFLOAT(M) B = A + 0.5D0 @@ -216,4 +218,3 @@ FUNCTION FMCHS(M,X) 999 FORMAT (24H0FMCHS DOES NOT CONVERGE,I6,1PE16.8) END C ----------------------- - diff --git a/lioamber/init_lio.f90 b/lioamber/init_lio.f90 index 4116e7275..8d444adc4 100644 --- a/lioamber/init_lio.f90 +++ b/lioamber/init_lio.f90 @@ -118,24 +118,24 @@ end subroutine lio_defaults !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! subroutine init_lio_common(natomin, Izin, nclatom, callfrom) - use garcha_mod, only : nunp, RMM, d, c, a, Nuc, ncont, cx, & - ax, Nucx, ncontx, cd, ad, Nucd, ncontd, indexii, & - indexiid, r, v, rqm, Em, Rm, pc, nnat, af, B, Iz, & - natom, ng0, ngd0, ngrid, nl, norbit, ntatom, & + use garcha_mod, only : nunp, RMM, d, c, a, r, v, rqm, Em, Rm, pc, nnat, B, & + Iz, natom, ng0, ngd0, ngrid, norbit, ntatom, & free_global_memory, little_cube_size,& assign_all_functions, energy_all_iterations, & remove_zero_weights, min_points_per_cube, & max_function_exponent, sphere_radius, M,Fock_Hcore, & Fock_Overlap, P_density, OPEN, timers, MO_coef_at, & - MO_coef_at_b, charge + MO_coef_at_b, charge, basis_set, fitting_set, Md, & + nl, basis, int_basis use ECP_mod, only : Cnorm, ecpmode use field_data, only : chrg_sq use fileio , only : lio_logo use fileio_data, only: style, verbose + use basis_subs, only: basis_init implicit none integer , intent(in) :: nclatom, natomin, Izin(natomin), callfrom - integer :: ng2, ngdDyn, ngDyn, ierr, ios, MM + integer :: ng2, ngdDyn, ngDyn, ierr, ios, MM, iostat if (verbose .gt. 2) then write(*,*) @@ -162,24 +162,18 @@ subroutine init_lio_common(natomin, Izin, nclatom, callfrom) ! Ngrid : n° of grid points (LS-SCF part). ! ! NOTES: Ngrid may be set to 0 in the case of Numerical Integration. For ! ! large systems, ng2 may result in <0 due to overflow. ! + if (.not. int_basis) basis_set = basis + call basis_init(basis_set, fitting_set, natom, Iz, iostat) + if (iostat .gt. 0) then + return + endif - ! Sets the dimensions for important arrays. - call DIMdrive(ngDyn,ngdDyn) - + ngDyn = M; ngdDyn = Md ng2 = 5*ngDyn*(ngDyn+1)/2 + 3*ngdDyn*(ngdDyn+1)/2 + & ngDyn + ngDyn*norbit + Ngrid - allocate(RMM(ng2) , d(natom, natom), c(ngDyn,nl) , a(ngDyn,nl) ,& - Nuc(ngDyn) , ncont(ngDyn) , cd(ngdDyn,nl) , ad(ngdDyn,nl) ,& - Nucd(ngdDyn), ncontd(ngdDyn) , indexii(ngDyn), indexiid(ngdDyn),& - v(ntatom,3) , Em(ntatom) , Rm(ntatom) , af(ngdDyn) ,& - nnat(200) , B(ngdDyn,3)) - - if (ngdDyn .gt. ngDyn) Then - allocate(cx(ngdDyn,nl), ax(ngdDyn,nl), Nucx(ngdDyn), ncontx(ngdDyn)) - else - allocate(cx(ngDyn,nl), ax(ngDyn,nl), Nucx(ngDyn), ncontx(ngDyn)) - endif + allocate(RMM(ng2) , d(natom, natom), v(ntatom,3), Em(ntatom), Rm(ntatom), & + nnat(200), B(ngdDyn,3)) ! Cnorm contains normalized coefficients of basis functions. ! Differentiate C for x^2,y^2,z^2 and xy,xz,yx (3^0.5 factor) @@ -336,15 +330,11 @@ subroutine init_lioamber_ehren(natomin, Izin, nclatom, charge_i, basis_i & , Fz_i, NBCH_i, propagator_i, writedens_i, tdrestart_i, dt_i & ) - use garcha_mod, only: M, first_step, doing_ehrenfest & - &, nshell, nuc, ncont, a, c, charge - use td_data, only: timedep, tdstep - - use basis_subs, only: basis_data_set, basis_data_norm - + use garcha_mod , only: M, first_step, doing_ehrenfest + use basis_subs , only: basis_setup_ehren + use td_data , only: timedep, tdstep use lionml_data, only: ndyn_steps, edyn_steps - - use liosubs, only: catch_error + use liosubs , only: catch_error implicit none @@ -373,13 +363,12 @@ subroutine init_lioamber_ehren(natomin, Izin, nclatom, charge_i, basis_i & , ntdstep_i, field_i, exter_i, a0_i, epsilon_i, Fx_i, Fy_i & , Fz_i, NBCH_i, propagator_i, writedens_i, tdrestart_i & ) - + call basis_setup_ehren() first_step=.true. if ( (ndyn_steps>0) .and. (edyn_steps>0) ) doing_ehrenfest=.true. - call basis_data_set(nshell(0),nshell(1),nshell(2),nuc,ncont,a,c) ! call basis_data_norm( M, size(c,2), c ) tdstep = (dt_i) * (41341.3733366d0) diff --git a/lioamber/lio_finalize.f b/lioamber/lio_finalize.f index 65b69670b..a4c399d7a 100644 --- a/lioamber/lio_finalize.f +++ b/lioamber/lio_finalize.f @@ -5,9 +5,14 @@ subroutine lio_finalize() ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! use garcha_mod - use ECP_mod, only : ecpmode - use fileio_data, only : style + use ECP_mod , only: ecpmode + use fileio_data, only: style + use basis_subs , only: basis_deinit + implicit none + call basis_deinit() + + ! This should not be here. if (dipole) then if (style) write(69,8703) CLOSE(69) @@ -24,12 +29,7 @@ subroutine lio_finalize() !--------------------------------------------------------------------! deallocate(r,v,rqm, Em, Rm) - deallocate(pc, Iz, cx, ax, cd, ad, c, a) - deallocate(Nuc,ncont,Nucx,ncontx,Nucd - > ,ncontd, indexii, indexiid, RMM, X) - deallocate(nnat, B, af) - deallocate(natomc,nnps,nnpp,nnpd,nns) - deallocate(nnd,nnp,atmin,jatc,d) + deallocate(pc, Iz, RMM, X,d) call g2g_timer_summary() call g2g_deinit() diff --git a/lioamber/liomods/garcha_mod.f b/lioamber/liomods/garcha_mod.f index 8f158d7ed..6d8dca169 100644 --- a/lioamber/liomods/garcha_mod.f +++ b/lioamber/liomods/garcha_mod.f @@ -10,7 +10,7 @@ module garcha_mod integer restart_freq, energy_freq real*8 GOLD, TOLD character*20 basis,whatis - character*40 basis_set, fitting_set + character*100 basis_set, fitting_set logical int_basis character*20 fcoord,fmulliken,frestart,frestartin,solv,solv2 logical exists,MEMO,predcoef diff --git a/liosolo/liomd.f90 b/liosolo/liomd.f90 index bbb39be54..fe836d9dc 100644 --- a/liosolo/liomd.f90 +++ b/liosolo/liomd.f90 @@ -8,7 +8,6 @@ program liomd use liosubs , only : write_energy, write_geom, set_masses, nuclear_verlet use fileio_data, only: verbose use ehrensubs - use basis_data implicit none real*8 :: escf, dipxyz(3), timeStep, Kenergy, Uenergy diff --git a/test/0_aguaCPU/basis b/test/0_aguaCPU/basis index c90dc3344..797cbbf43 100644 --- a/test/0_aguaCPU/basis +++ b/test/0_aguaCPU/basis @@ -18,8 +18,8 @@ gaussian 0.1938135 1.0000000 0.8000000 1.0000000 8 13 13 - 1 1 1 1 1 1 1 1 1 1 1 1 1 - 0 0 0 0 0 0 0 1 1 1 2 2 2 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 0 0 0 0 0 0 0 1 1 1 2 2 2 2000.00000000 1.00000000 400.00000000 1.00000000 100.00000000 1.00000000 diff --git a/test/0_aguaGPU/basis b/test/0_aguaGPU/basis index c90dc3344..312e42bca 100644 --- a/test/0_aguaGPU/basis +++ b/test/0_aguaGPU/basis @@ -18,8 +18,8 @@ gaussian 0.1938135 1.0000000 0.8000000 1.0000000 8 13 13 - 1 1 1 1 1 1 1 1 1 1 1 1 1 - 0 0 0 0 0 0 0 1 1 1 2 2 2 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 0 0 0 0 0 0 0 1 1 1 2 2 2 2000.00000000 1.00000000 400.00000000 1.00000000 100.00000000 1.00000000 diff --git a/test/1_Oxy-molGPU/DZVP b/test/1_Oxy-molGPU/DZVP index 4d54c14e2..b97bfa3ca 100644 --- a/test/1_Oxy-molGPU/DZVP +++ b/test/1_Oxy-molGPU/DZVP @@ -83,8 +83,8 @@ gaussian 0.5115279 -0.4862694 0.0900000 1.0000000 26 20 20 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 0.44000000E+05 1.00000000 0.08800000E+05 1.00000000 0.22000000E+04 1.00000000 @@ -107,8 +107,8 @@ gaussian 0.13500000E+00 1.00000000 gaussian 7 15 6 -6 2 1 4 1 1 -0 0 0 1 1 2 + 6 2 1 4 1 1 + 0 0 0 1 1 2 3845.4149000 0.0020186 577.5332300 0.0154078 131.3198300 0.0753714 @@ -142,8 +142,8 @@ gaussian 0.32000000 1.00000000 gaussian 16 21 8 -6 3 2 1 5 2 1 1 -0 0 0 0 1 1 1 2 + 6 3 2 1 5 2 1 1 + 0 0 0 0 1 1 1 2 23050.0670000 -0.0017567 3451.8663000 -0.0134894 781.2786700 -0.0674325 @@ -166,8 +166,8 @@ gaussian 0.1021157 1.0000000 0.6500000 1.0000000 16 17 17 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 0 0 0 0 0 0 0 0 0 1 1 1 1 2 2 2 2 0.16384000E+05 1.00000000 0.32770000E+04 1.00000000 0.81900000E+03 1.00000000 @@ -205,8 +205,8 @@ gaussian 0.1938135 1.0000000 0.8000000 1.0000000 8 13 13 - 1 1 1 1 1 1 1 1 1 1 1 1 1 - 0 0 0 0 0 0 0 1 1 1 2 2 2 + 1 1 1 1 1 1 1 1 1 1 1 1 1 + 0 0 0 0 0 0 0 1 1 1 2 2 2 2000.00000000 1.00000000 400.00000000 1.00000000 100.00000000 1.00000000