comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 22:25:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Correct spin unrestricted handling in xTB [ci skip]

Browse source
This commit is contained in:
Daniel Mejia-Rodriguez 2022-04-08 14:24:10 -04:00
parent 9e219fe8cc
commit e1bb7abb30
10 changed files with 330 additions and 207 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/xtb/GNUmakefile
View file

@ -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

8
src/xtb/xtb.fh
View file

@ -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

160
src/xtb/xtb_energy.F
View file

@ -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

163
src/xtb/xtb_energy_gradient.F
View file

@ -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

2
src/xtb/xtb_guess.F
View file

@ -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)

34
src/xtb/xtb_input.F
View file

@ -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

69
src/xtb/xtb_print_param.F Normal file
View file

@ -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

7
src/xtb/xtb_rdinput.F
View file

@ -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

42
src/xtb/xtb_spinpol.F Normal file
View file

@ -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

48
src/xtb/xtb_utils.F
View file

@ -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