diff --git a/src/xtb/GNUmakefile b/src/xtb/GNUmakefile index a796887c1c..be8510584d 100644 --- a/src/xtb/GNUmakefile +++ b/src/xtb/GNUmakefile @@ -2,8 +2,8 @@ ifndef USE_TBLITE $(error USE_TBLITE must be set in order to use xTB methods) endif -OBJ = xtb_energy.o xtb_input.o xtb_utils.o xtb_energy_gradient.o \ - xtb_rdinput.o xtb_guess.o +OBJ = xtb_energy.o xtb_input.o xtb_energy_gradient.o \ + xtb_rdinput.o xtb_guess.o xtb_spinpol.o xtb_print_param.o LIBRARY = libnwxtb.a diff --git a/src/xtb/xtb.fh b/src/xtb/xtb.fh index 6d1ee32518..15635ca2c3 100644 --- a/src/xtb/xtb.fh +++ b/src/xtb/xtb.fh @@ -1,7 +1,7 @@ - double precision acc, rcharge, kt + double precision acc, rcharge, kt, broyden logical gfn1, gfn2, gff, guess - integer uhf, verb + integer uhf, verb, nspin common /cxtb/ - $ acc, rcharge, kt, + $ acc, rcharge, kt, broyden, $ gfn1, gfn2, gff, guess, - $ uhf, verb + $ uhf, verb, nspin diff --git a/src/xtb/xtb_energy.F b/src/xtb/xtb_energy.F index cf045dae3e..cdf54651d1 100644 --- a/src/xtb/xtb_energy.F +++ b/src/xtb/xtb_energy.F @@ -14,7 +14,6 @@ use tblite_xtb_gfn2, only : new_gfn2_calculator, export_gfn2_param use tblite_xtb_gfn1, only : new_gfn1_calculator, export_gfn1_param use tblite_xtb_singlepoint, only : xtb_singlepoint - use tblite_scf_info, only : scf_info,atom_resolved,shell_resolved use tblite_param, only : param_record #ifdef USE_OPENMP use omp_lib, only : omp_get_max_threads @@ -46,22 +45,20 @@ character(len=8) theory character(len=255) temp - integer*4 :: iuhf4, one4, zero4, nspin4, verb4 - integer*4,allocatable :: atnum(:) + integer*4 :: iuhf4, nspin4, verb4 + integer*4,allocatable :: atnum4(:) double precision energy double precision, allocatable :: coords(:,:) double precision, allocatable :: rcov(:), trans(:,:), cn(:) - type(structure_type) :: mol - type(xtb_calculator) :: calc - type(wavefunction_type) :: wfn - type(context_type) :: ctx - type(param_record) :: param - type(scf_info) :: info - - parameter(one4=1, zero4=0) + type(structure_type),allocatable :: mol + type(xtb_calculator),allocatable :: calc + type(wavefunction_type),allocatable :: wfn + type(context_type),allocatable :: ctx + type(param_record),allocatable :: param + #ifdef USE_OPENMP integer maxthreads maxthreads = omp_get_max_threads() @@ -70,10 +67,11 @@ xtb_energy = .false. status = rtdb_parallel(.true.) ! Broadcast reads to all processes + + call util_print_push() call util_print_rtdb_load(rtdb,'xtb') - oprint = util_print('information', print_low) .and. - $ ga_nodeid().eq.0 + oprint = util_print('information',print_low).and.ga_nodeid().eq.0 oparam = ga_nodeid().eq.0 .and.util_print('params',print_high) ovecs = ga_nodeid().eq.0 .and. util_print('vecs',print_default) @@ -91,7 +89,7 @@ ! ! Initialize ! - if (.not.rtdb_put(rtdb, 'xtb:energy', mt_dbl, 1, 0d0)) + if (.not.rtdb_put(rtdb, 'xtb:energy', mt_dbl, 1, 0.0d0)) $ call errquit('xtb_energy: faild to store energy',0,RTDB_ERR) if (.not.rtdb_put(rtdb, 'xtb:converged', mt_log, 1, .false.)) $ call errquit('xtb_energy: failed to store converged',0,RTDB_ERR) @@ -110,8 +108,7 @@ verb4 = verb iuhf4 = uhf - nspin4 = 1 -C if (iuhf4.gt.0) nspin4 = 2 + nspin4 = nspin ! ! Allocate space for atomic coordinates @@ -156,22 +153,42 @@ C if (iuhf4.gt.0) nspin4 = 2 if (nat.ne.ncenters) $ call errquit('xtb_energy: not ready for dummy atoms',0,GEOM_ERR) - allocate(atnum(ncenters),stat=istat) + allocate(atnum4(ncenters),stat=istat) if (istat.ne.0) $ call errquit('xtb_energy: could not allocate numbers',0,MA_ERR) do i=0,ncenters-1 - atnum(i+1) = nint(dbl_mb(icharge+i)) + atnum4(i+1) = nint(dbl_mb(icharge+i)) enddo energy = 0.0d0 + ! + ! Line needed for MOLDEN + ! if (oprint) then write(luout,'(10X,A,I6)') 'No. of atoms :',ncenters write(luout,*) endif + + ! + ! Allocate xTB objects + ! + allocate(mol) + allocate(calc) + allocate(wfn) + allocate(ctx) + allocate(param) + if (ga_nodeid().eq.0) then - call new(mol,atnum,coords,charge=rcharge,uhf=iuhf4) + ! + ! Pass geometry info and create mol object + ! + call new(mol,atnum4,coords,charge=rcharge,uhf=iuhf4) + + ! + ! Obtain the appropriate calculator + ! if (gfn2) then call new_gfn2_calculator(calc, mol) call export_gfn2_param(param) @@ -179,72 +196,64 @@ C if (iuhf4.gt.0) nspin4 = 2 call new_gfn1_calculator(calc, mol) call export_gfn1_param(param) endif + + ! + ! Set other parameters + ! + calc%mixer_damping = broyden + + ! + ! Initialize wavefunction object + ! call new_wavefunction(wfn,mol%nat,calc%bas%nsh,calc%bas%nao, $ nspin4,kt) + ! + ! Setup spin-unrestricted calculations + ! + if (nspin.eq.2) call xtb_spinpol(mol, calc, calc%bas) + ! ! Compute starting charges + ! if (guess) then call xtb_guess(mol,wfn%qat,wfn%qsh,param,mol%nat, - $ calc%bas%nsh,calc%bas%nsh_at,atnum) + $ calc%bas%nsh,calc%bas%nsh_at,atnum4) + ! + ! Copy charges to second spin-channel, if available + ! + if (nspin.eq.2) then + wfn%qat(:,2) = wfn%qat(:,1) + wfn%qsh(:,2) = wfn%qsh(:,1) + endif + endif + ! + ! Print all parameters + ! if (oparam) then - - info = calc%coulomb%es3%variable_info() - select case(info%charge) - case (atom_resolved) - temp = 'atom-resolved' - case (shell_resolved) - temp = 'shell-resolved' - end select - - - 5000 format(/,10X,'No. of atomic orbitals: ',I10, - $ /,10X,'No. of shells : ',I10, - $ /,10X,'Electronic temp. : ',F10.6, - $ /,10X,'Accuracy : ',F10.6, - $ /,10X,' Integral cutoff : ',F10.6, - $ /,10X,' E. convergence : ',F10.6, - $ /,10X,' W.F. convergence : ',F10.6, - $ /,10X,'Broyden damping : ',F10.6, - $ /,10X,'Couomb 3rd order : ',A, - $ /,10X,'H0-scaling : ',3(F10.6), - $ /,10X,'s8 : ',F10.6, - $ /,10X,'a1 : ',F10.6, - $ /,10X,'a2 : ',F10.6, - $ /,10X,'s9 : ',F10.6, - $ /,10X,'kExp : ',2(F10.6), - $ /,10X,'alpha : ',F10.6, - $ /,10X,'a3 : ',F10.6, - $ /,10X,'a5 : ',F10.6, - $ /,10X,'cn-shift : ',F10.6, - $ /,10X,'cn-exp : ',F10.6, - $ /,10X,'max-rad : ',F10.6) - - - write(luout,5000) calc%bas%nao,calc%bas%nsh,kt/kb_au,acc, - $ calc%bas%intcut, 1d-6*acc, 2d-5*acc,calc%mixer_damping, - $ temp(1:len(trim(temp))), - $ (param%hamiltonian%ksh(i,i),i=0,param%hamiltonian%lmax), - $ param%dispersion%s8,param%dispersion%a1,param%dispersion%a2, - $ param%dispersion%s9,param%repulsion%kexp, - $ param%repulsion%klight,param%charge%gexp, - $ param%multipole%dmp3,param%multipole%dmp5, - $ param%multipole%shift,param%multipole%kexp, - $ param%multipole%rmax - call util_flush(luout) + call xtb_print_param(calc, param) endif + ! + ! Actually do the sinlgepoint calculations + ! call xtb_singlepoint(ctx,mol,calc,wfn,acc,energy, $ verbosity=verb4) + ! + ! Total Energy Output + ! write(luout,*) write(luout,'(A,F15.8)') ' Total xTB energy = ',energy write(luout,*) + ! + ! Print orbital energies + ! if (ovecs) then write(luout,*) write(luout,*) ' Alpha Orbital Energies ' @@ -266,17 +275,29 @@ C if (iuhf4.gt.0) nspin4 = 2 write(luout,*) endif endif - - call xtb_deallocate_wavefunction(wfn) - call xtb_deallocate_structure(mol) endif - + ! + ! Check if something failed during the calculation + ! + xtb_energy = .not.ctx%failed() + call ga_brdcst(101,xtb_energy,ma_sizeof(mt_log,1,mt_byte),0) + + ! + ! Deallocate xTB objects + ! + deallocate(ctx) + deallocate(wfn) + deallocate(calc) + deallocate(mol) + deallocate(param) + + if (.not.ma_chop_stack(lcharge)) $ call errquit('xtb_energy: failed to pop stack',0,MA_ERR) if (.not.geom_destroy(geom)) $ call errquit('xtb_energy: could not destroy geom',0,GEOM_ERR) - deallocate(coords,atnum,stat=istat) + deallocate(coords,atnum4,stat=istat) if (istat.ne.0) $ call errquit('xtb_energy: failed to deallocate',0,MA_ERR) @@ -286,7 +307,6 @@ C if (iuhf4.gt.0) nspin4 = 2 if (.not.rtdb_put(rtdb,'xtb:converged',mt_log,1,.true.)) $ call errquit('xtb_energy: failed to put in rtdb',0,RTDB_ERR) - xtb_energy = .not.ctx%failed() call util_print_pop() #ifdef USE_OPENMP diff --git a/src/xtb/xtb_energy_gradient.F b/src/xtb/xtb_energy_gradient.F index 556a77652e..23b0b0f624 100644 --- a/src/xtb/xtb_energy_gradient.F +++ b/src/xtb/xtb_energy_gradient.F @@ -14,7 +14,6 @@ use tblite_xtb_gfn2, only : new_gfn2_calculator, export_gfn2_param use tblite_xtb_gfn1, only : new_gfn1_calculator, export_gfn1_param use tblite_xtb_singlepoint, only : xtb_singlepoint - use tblite_scf_info, only : scf_info,atom_resolved,shell_resolved use tblite_param, only : param_record #ifdef USE_OPENMP use omp_lib, only : omp_get_max_threads @@ -50,20 +49,18 @@ character(len=8) theory character(len=255) temp - integer*4 :: iuhf4, one4, zero4, nspin4 - integer*4,allocatable :: atnum(:) + integer*4 :: iuhf4, verb4, nspin4 + integer*4,allocatable :: atnum4(:) double precision energy, sigma(3,3), dipole(3) double precision, allocatable :: coords(:,:),grad(:,:) double precision, allocatable :: rcov(:), trans(:,:), cn(:) - type(structure_type) :: mol - type(xtb_calculator) :: calc - type(wavefunction_type) :: wfn - type(context_type) :: ctx - type(scf_info) :: info - type(param_record) :: param - parameter(one4=1, zero4=0) + type(structure_type),allocatable :: mol + type(xtb_calculator),allocatable :: calc + type(wavefunction_type),allocatable :: wfn + type(context_type),allocatable :: ctx + type(param_record),allocatable :: param #ifdef USE_OPENMP integer maxthreads @@ -73,11 +70,13 @@ xtb_energy_gradient = .false. status = rtdb_parallel(.true.) ! Broadcast reads to all processes + + ! + ! Printing options + ! call util_print_push() call util_print_rtdb_load(rtdb,'xtb') - - oprint = util_print('information', print_low) .and. - $ ga_nodeid().eq.0 + oprint = util_print('information',print_low).and.ga_nodeid().eq.0 oparam = ga_nodeid().eq.0 .and.util_print('params',print_high) ovecs = ga_nodeid().eq.0 .and. util_print('vecs',print_default) @@ -94,7 +93,7 @@ ! ! Initialize ! - if (.not.rtdb_put(rtdb, 'xtb:energy', mt_dbl, 1, 0d0)) + if (.not.rtdb_put(rtdb, 'xtb:energy', mt_dbl, 1, 0.0d0)) $ call errquit('xtb_energy: faild to store energy',0,RTDB_ERR) if (.not.rtdb_put(rtdb, 'xtb:converged', mt_log, 1, .false.)) $ call errquit('xtb_energy: failed to store converged',0,RTDB_ERR) @@ -110,10 +109,11 @@ $ call errquit('xtb_energy: geom_ncent failed',0,GEOM_ERR) if (.not.rtdb_get(rtdb,'charge',mt_dbl,1,rcharge)) $ rcharge = 0d0 - iuhf4 = uhf - nspin4 = 1 -C if (iuhf4.gt.0) nspin4 = 2 + + iuhf4 = uhf + nspin4 = nspin + verb4 = verb C symmetry or not if (.not. rtdb_get(rtdb, 'gradients:use symmetry', mt_log, 1, @@ -178,30 +178,48 @@ c if (nat.ne.ncenters) $ call errquit('xtb_energy: not ready for dummy atoms',0,GEOM_ERR) - allocate(atnum(ncenters),stat=istat) + ! + ! Atomic numbers + ! + allocate(atnum4(ncenters),stat=istat) if (istat.ne.0) $ call errquit('xtb_energy: could not allocate numbers',0,MA_ERR) - - ! Atomic numbers do i=0,ncenters-1 - atnum(i+1) = nint(dbl_mb(icharge+i)) + atnum4(i+1) = nint(dbl_mb(icharge+i)) enddo energy = 0.0d0 grad(:,:) = 0.0d0 sigma(:,:) = 0.0d0 + ! + ! This line is needed for MOLDEN + ! if (oprint) then write(luout,'(10X,A,I6)') 'No. of atoms :',ncenters write(luout,*) endif + ! + ! Allocate xTB objects + ! + allocate(mol) + allocate(param) + allocate(wfn) + allocate(calc) + allocate(ctx) + if (ga_nodeid().eq.0) then + status = rtdb_parallel(.false.) + ! ! Create mol object - call new(mol,atnum,coords,charge=rcharge,uhf=iuhf4) + ! + call new(mol,atnum4,coords,charge=rcharge,uhf=iuhf4) + ! ! Select the calculator type and get params + ! if (gfn2) then theory = 'xTB GFN2' call new_gfn2_calculator(calc, mol) @@ -212,81 +230,65 @@ c call export_gfn1_param(param) endif + ! ! Initialize wavefunction object + ! call new_wavefunction(wfn,mol%nat,calc%bas%nsh,calc%bas%nao, $ nspin4,kt) + ! ! Compute starting charges + ! status = rtdb_get(rtdb,'xtb:qat',mt_dbl,nat,wfn%qat) status = status .and. rtdb_get(rtdb,'xtb:qsh',mt_dbl, $ int(calc%bas%nsh),wfn%qsh) if (guess.and..not.status) then call xtb_guess(mol,wfn%qat,wfn%qsh,param,mol%nat, - $ calc%bas%nsh,calc%bas%nsh_at,atnum) + $ calc%bas%nsh,calc%bas%nsh_at,atnum4) endif + ! + ! Setup spin-unrestricted calculation + ! + if (nspin.eq.2) then + call xtb_spinpol(mol, calc, calc%bas) + endif + ! + ! Set other parameters + ! + calc%mixer_damping = broyden + + ! + ! Print params + ! if (oparam) then - - info = calc%coulomb%es3%variable_info() - select case(info%charge) - case (atom_resolved) - temp = 'atom-resolved' - case (shell_resolved) - temp = 'shell-resolved' - end select - - - 5000 format(/,10X,'No. of atomic orbitals: ',I10, - $ /,10X,'No. of shells : ',I10, - $ /,10X,'Electronic temp. : ',F10.6, - $ /,10X,'Accuracy : ',F10.6, - $ /,10X,' Integral cutoff : ',F10.6, - $ /,10X,' E. convergence : ',F10.6, - $ /,10X,' W.F. convergence : ',F10.6, - $ /,10X,'Broyden damping : ',F10.6, - $ /,10X,'Couomb 3rd order : ',A, - $ /,10X,'H0-scaling : ',3(F10.6), - $ /,10X,'s8 : ',F10.6, - $ /,10X,'a1 : ',F10.6, - $ /,10X,'a2 : ',F10.6, - $ /,10X,'s9 : ',F10.6, - $ /,10X,'kExp : ',2(F10.6), - $ /,10X,'alpha : ',F10.6, - $ /,10X,'a3 : ',F10.6, - $ /,10X,'a5 : ',F10.6, - $ /,10X,'cn-shift : ',F10.6, - $ /,10X,'cn-exp : ',F10.6, - $ /,10X,'max-rad : ',F10.6) - - - write(luout,5000) calc%bas%nao,calc%bas%nsh,kt/kb_au,acc, - $ calc%bas%intcut, 1d-6*acc, 2d-5*acc,calc%mixer_damping, - $ temp(1:len(trim(temp))), - $ (param%hamiltonian%ksh(i,i),i=0,param%hamiltonian%lmax), - $ param%dispersion%s8,param%dispersion%a1,param%dispersion%a2, - $ param%dispersion%s9,param%repulsion%kexp, - $ param%repulsion%klight,param%charge%gexp, - $ param%multipole%dmp3,param%multipole%dmp5, - $ param%multipole%shift,param%multipole%kexp, - $ param%multipole%rmax - call util_flush(luout) + call xtb_print_param(calc, param) endif + ! + ! Obtain energy and gradients + ! call xtb_singlepoint(ctx,mol,calc,wfn,acc,energy,grad,sigma, - $ verbosity=int(verb,kind=4)) + $ verbosity=verb4) + ! + ! Total energy output + ! write(luout,*) write(luout,'(A,F15.8)') ' Total xTB energy = ',energy write(luout,*) -c -C zero force contributions on inactive atoms + ! + ! zero force contributions on inactive atoms + ! call grad_active_atoms(rtdb, nat, log_mb(iact), nactive) call zero_forces(grad, log_mb(iact), nat) - -C symmetrize + + ! + ! symmetrize + ! if (oskel) call sym_grad_symmetrize(geom, grad) ! @@ -343,17 +345,25 @@ C symmetrize dipole(1:3) = sum(wfn%dpat(:,:,1),dim=2) status = rtdb_put(rtdb,'task:dipole',mt_dbl,3,dipole) - - call xtb_deallocate_wavefunction(wfn) - call xtb_deallocate_structure(mol) + status = rtdb_parallel(.true.) endif + xtb_energy_gradient = .not.ctx%failed() + + ! + ! Deallocate all xTB objects + ! + deallocate(ctx) + deallocate(wfn) + deallocate(calc) + deallocate(mol) + deallocate(param) if (.not.ma_chop_stack(lcharge)) $ call errquit('xtb_energy: failed to pop stack',0,MA_ERR) if (.not.geom_destroy(geom)) $ call errquit('xtb_energy: could not destroy geom',0,GEOM_ERR) - deallocate(coords,atnum,stat=istat) + deallocate(coords,atnum4,stat=istat) if (istat.ne.0) $ call errquit('xtb_energy: failed to deallocate',0,MA_ERR) @@ -374,7 +384,6 @@ C symmetrize $ call errquit('xtb_energy: failed to deallocate',0,MA_ERR) - xtb_energy_gradient = .not.ctx%failed() call util_print_pop() #ifdef USE_OPENMP diff --git a/src/xtb/xtb_guess.F b/src/xtb/xtb_guess.F index 5f9c2ff543..bb1c71d2b6 100644 --- a/src/xtb/xtb_guess.F +++ b/src/xtb/xtb_guess.F @@ -30,7 +30,7 @@ call model%solve(mol, cn, qvec=qat) ! transform atomic charges into shell charges - ksh = 0 + ksh = 0 do i=1,nat tot = sum(param%record(atnum(i))%refocc(:)) do j=1,nshl_at(i) diff --git a/src/xtb/xtb_input.F b/src/xtb/xtb_input.F index 5446d03480..1ad3eaacd5 100644 --- a/src/xtb/xtb_input.F +++ b/src/xtb/xtb_input.F @@ -24,13 +24,14 @@ c integer rtdb character*255 test c - integer, parameter :: num_dirs = 9 + integer, parameter :: num_dirs = 11 character*11 dirs(num_dirs) data dirs /'acc','uhf', 'method', 'etemp', 'guess', - $ 'verbosity', 'print', 'noprint', 'end'/ + $ 'verbosity', 'print', 'noprint', 'nspin', + $ 'broyden', 'end'/ c - integer uhf, ind, verb - double precision acc, etemp + integer uhf, ind, verb, nspin + double precision acc, etemp, broyden character(len=4) method c @@ -48,7 +49,7 @@ c call errquit('xtb_input: unknown directive', 0, INPUT_ERR) end if c - goto (100, 200, 300, 400, 500, 600, 700, 800, 9999) ind + goto (100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 9999) ind call errquit('xtb_input: unknown directive',ind, INPUT_ERR) c @@ -116,6 +117,29 @@ c 800 call util_print_input(rtdb, 'xtb') goto 10 c +c spin-polarization treatment +c + 900 if(.not.inp_i(nspin)) + $ call errquit('xtb_input: nspin value not found',900,INPUT_ERR) + if (.not.(nspin.eq.1 .or. nspin.eq.2)) then + call errquit('xtb_input: nspin can be 1 or 2',900,INPUT_ERR) + endif + if (.not.rtdb_put(rtdb, 'xtb:nspin', mt_int, 1, nspin)) + $ call errquit('xtb_input: rtdb_put failed',900,RTDB_ERR) + goto 10 +c +c broyden mixer parameter +c + 1000 if(.not.inp_f(broyden)) + $ call errquit('xtb_input: broyden mixer value not found',1000, + $ INPUT_ERR) + if (broyden.le.0.0d0 .or. broyden .gt. 1.0d0) + $ call errquit('xtb_input: broyden value should lie (0,1]',1001, + $ INPUT_ERR) + if (.not.rtdb_put(rtdb, 'xtb:broyden', mt_dbl, 1, broyden)) + $ call errquit('xtb_input: rtdb_put failed',1002,RTDB_ERR) + goto 10 +c c normal termination c 9999 return diff --git a/src/xtb/xtb_print_param.F b/src/xtb/xtb_print_param.F new file mode 100644 index 0000000000..5b17cba2d3 --- /dev/null +++ b/src/xtb/xtb_print_param.F @@ -0,0 +1,69 @@ + subroutine xtb_print_param(calc, param) + use tblite_xtb_calculator, only: xtb_calculator + use tblite_scf_info, only : scf_info,atom_resolved,shell_resolved + use tblite_param, only : param_record + + implicit none + +#include "stdio.fh" +#include "util.fh" +#include "util_params.fh" +#include "xtb.fh" + + type(xtb_calculator) :: calc + type(param_record) :: param + type(scf_info),allocatable :: info + + character(len=255) temp + + integer i + + allocate(info) + info = calc%coulomb%es3%variable_info() + + select case(info%charge) + case (atom_resolved) + temp = 'atom-resolved' + case (shell_resolved) + temp = 'shell-resolved' + end select + + + 5000 format(/,10X,'No. of atomic orbitals: ',I10, + $ /,10X,'No. of shells : ',I10, + $ /,10X,'Electronic temp. : ',F10.6, + $ /,10X,'Accuracy : ',F10.6, + $ /,10X,' Integral cutoff : ',F10.6, + $ /,10X,' E. convergence : ',F10.6, + $ /,10X,' W.F. convergence : ',F10.6, + $ /,10X,'Broyden damping : ',F10.6, + $ /,10X,'Couomb 3rd order : ',A, + $ /,10X,'H0-scaling : ',3(F10.6), + $ /,10X,'s8 : ',F10.6, + $ /,10X,'a1 : ',F10.6, + $ /,10X,'a2 : ',F10.6, + $ /,10X,'s9 : ',F10.6, + $ /,10X,'kExp : ',2(F10.6), + $ /,10X,'alpha : ',F10.6, + $ /,10X,'a3 : ',F10.6, + $ /,10X,'a5 : ',F10.6, + $ /,10X,'cn-shift : ',F10.6, + $ /,10X,'cn-exp : ',F10.6, + $ /,10X,'max-rad : ',F10.6) + + + write(luout,5000) calc%bas%nao,calc%bas%nsh,kt/kb_au,acc, + $ calc%bas%intcut, 1d-6*acc, 2d-5*acc,calc%mixer_damping, + $ temp(1:len(trim(temp))), + $ (param%hamiltonian%ksh(i,i),i=0,param%hamiltonian%lmax), + $ param%dispersion%s8,param%dispersion%a1,param%dispersion%a2, + $ param%dispersion%s9,param%repulsion%kexp, + $ param%repulsion%klight,param%charge%gexp, + $ param%multipole%dmp3,param%multipole%dmp5, + $ param%multipole%shift,param%multipole%kexp, + $ param%multipole%rmax + + call util_flush(luout) + deallocate(info) + + end diff --git a/src/xtb/xtb_rdinput.F b/src/xtb/xtb_rdinput.F index c86a705c84..4a9ec21570 100644 --- a/src/xtb/xtb_rdinput.F +++ b/src/xtb/xtb_rdinput.F @@ -16,7 +16,10 @@ if (.not.rtdb_get(rtdb,'xtb:uhf',mt_int,1,uhf)) uhf=0 if (.not.rtdb_cget(rtdb,'xtb:method',1,method)) method='gfn2' if (.not.rtdb_get(rtdb,'xtb:verb',mt_int,1,verb)) verb=1 + if (.not.rtdb_get(rtdb,'xtb:nspin',mt_int,1,nspin)) nspin=1 if (.not.rtdb_get(rtdb,'xtb:guess',mt_log,1,guess)) guess=.false. + if (.not.rtdb_get(rtdb,'xtb:broyden',mt_dbl,1,broyden)) + $ broyden=0.4d0 if (method.eq.'gfn2') then gfn2 = .true. @@ -44,5 +47,9 @@ $ call errquit('xtb_rdinput: cannot put in rtdb',5,RTDB_ERR) if (.not.rtdb_cput(rtdb,'xtb:method',1,method)) $ call errquit('xtb_rdinput: cannot put in rtdb',6,RTDB_ERR) + if (.not.rtdb_put(rtdb,'xtb:nspin',mt_int,1,nspin)) + $ call errquit('xtb_rdinput: cannot put in rtdb',7,RTDB_ERR) + if (.not.rtdb_put(rtdb,'xtb:broyden',mt_dbl,1,broyden)) + $ call errquit('xtb_rdinput: cannot put in rtdb',8,RTDB_ERR) end subroutine diff --git a/src/xtb/xtb_spinpol.F b/src/xtb/xtb_spinpol.F new file mode 100644 index 0000000000..25db1c0e61 --- /dev/null +++ b/src/xtb/xtb_spinpol.F @@ -0,0 +1,42 @@ + subroutine xtb_spinpol(mol,calc,bas) + use mctc_io, only: structure_type + use tblite_container, only : container_type + use tblite_spin, only : spin_polarization, new_spin_polarization + use tblite_data_spin, only : get_spin_constant + use tblite_basis_type, only : basis_type + use tblite_xtb_calculator, only: xtb_calculator + + implicit none + + type(structure_type) :: mol + type(xtb_calculator) :: calc + type(basis_type) :: bas + type(spin_polarization), allocatable :: spin + class(container_type), allocatable :: cont + + integer i, ish, jsh + double precision, allocatable :: wll(:,:,:) + + allocate(spin) + allocate(wll(bas%nsh,bas%nsh,mol%nid),source=0.0d0) + do i=1,mol%nid + do ish=1,bas%nsh_id(i) + do jsh=1,bas%nsh_id(i) + wll(ish,jsh,i) = get_spin_constant( + $ bas%cgto(ish,i)%ang, + $ bas%cgto(jsh,i)%ang, + $ mol%num(i)) + enddo + enddo + enddo + call new_spin_polarization(spin, mol, wll, bas%nsh_id) + call move_alloc(spin, cont) + call calc%push_back(cont) + + ! + ! Both xTB objects should be deallocated by now + ! + if (allocated(spin)) deallocate(spin) + if (allocated(cont)) deallocate(cont) + + end subroutine diff --git a/src/xtb/xtb_utils.F b/src/xtb/xtb_utils.F deleted file mode 100644 index 6271f4fc8c..0000000000 --- a/src/xtb/xtb_utils.F +++ /dev/null @@ -1,48 +0,0 @@ - subroutine xtb_deallocate_wavefunction(wfn) - use tblite_wavefunction_type, only : wavefunction_type - implicit none - type(wavefunction_type) :: wfn - - deallocate(wfn%homo) - deallocate(wfn%nel) - deallocate(wfn%n0at) - deallocate(wfn%n0sh) - deallocate(wfn%density) - deallocate(wfn%coeff) - deallocate(wfn%emo) - deallocate(wfn%focc) - deallocate(wfn%qat) - deallocate(wfn%qsh) - deallocate(wfn%dpat) - deallocate(wfn%qpat) - - end subroutine - - subroutine xtb_deallocate_structure(struct) - use mctc_io, only: structure_type - implicit none - type(structure_type) :: struct - - deallocate(struct%id) - deallocate(struct%xyz) - deallocate(struct%lattice) - deallocate(struct%periodic) - deallocate(struct%num) - deallocate(struct%sym) - - end - - subroutine xtb_deallocate_params(param) - use tblite_param, only : param_record - implicit none - type(param_record) :: param - deallocate(param%name) - deallocate(param%reference) - deallocate(param%dispersion) - deallocate(param%repulsion) - deallocate(param%charge) - deallocate(param%multipole) - deallocate(param%halogen) - deallocate(param%thirdorder) - deallocate(param%record) - end