mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 06:05:44 -04:00
routines cleanup
This commit is contained in:
Edoardo Apra
parent
28d1185904
commit
449dce3138
4 changed files with 568 additions and 485 deletions
|
|
@ -1,12 +1,12 @@
|
|||
#
|
||||
# $Id: GNUmakefile,v 1.48 2000-03-27 19:19:23 edo Exp $
|
||||
# $Id: GNUmakefile,v 1.49 2000-03-30 03:07:30 edo Exp $
|
||||
#
|
||||
|
||||
OBJ = dft_grad_top.o dftg_force.o \
|
||||
dftg_getxc.o dft_gradients.o \
|
||||
grad_dft.o dftg_quadv0.o \
|
||||
dftg_getcden.o dftg_getxden.o \
|
||||
dftg_bal_sum.o dftg_gridv0a.o dftg_gridv0.o
|
||||
dftg_bal_sum.o dftg_gridv0a.o dftg_gridv0.o dftg_gridv0b.o
|
||||
OBJ_OPTIMIZE = dftg_cdfit.o dftg_quadv0a.o dftg_quadv0b.o dftg_ssw.o\
|
||||
dftg_grad2.o dftg_hf_x.o dint_block_j.o j_nucder_gen.o \
|
||||
dftg_beckew.o dftg_wderiv.o
|
||||
|
|
@ -16,11 +16,11 @@ ifeq ($(NWCHEM_TARGET),SGITFP)
|
|||
grad_dft.o dftg_quadv0.o dftg_gridv0a.o\
|
||||
dftg_getcden.o dftg_getxden.o\
|
||||
dftg_bal_sum.o dftg_hf_x.o dftg_quadv0a.o dftg_beckew.o dftg_wderiv.o\
|
||||
dint_block_j.o j_nucder_gen.o dftg_ssw.o dftg_gridv0.o
|
||||
dint_block_j.o j_nucder_gen.o dftg_ssw.o dftg_gridv0.o dftg_gridv0b.o
|
||||
OBJ_OPTIMIZE = dftg_quadv0b.o
|
||||
endif
|
||||
|
||||
USES_BLAS = dftg_cdfit.F dftg_quadv0a.F dftg_quadv0b.F dftg_wderiv.F j_nucder_gen.F dft_gradients.F dftg_force.F dftg_gridv0a.F
|
||||
USES_BLAS = dftg_cdfit.F dftg_quadv0a.F dftg_quadv0b.F dftg_wderiv.F j_nucder_gen.F dft_gradients.F dftg_force.F
|
||||
|
||||
LIB_INCLUDES = -I../nwdft/include -I../nwdft/grid
|
||||
|
||||
|
|
|
|||
|
|
@ -4,24 +4,24 @@ c routine a wrapper for the more general routine
|
|||
c
|
||||
c BGJ - 8/98
|
||||
c
|
||||
Subroutine dftg_gridv0(rtdb, iga_dens, force, natom, oactive,
|
||||
Subroutine dftg_gridv0(rtdb, iga_dens, force, natoms, oactive,
|
||||
& nactive, rdens_atom, cntoce, cntobfr,
|
||||
& cetobfr)
|
||||
c
|
||||
c $Id: dftg_gridv0.F,v 1.1 2000-03-22 03:02:16 edo Exp $
|
||||
c $Id: dftg_gridv0.F,v 1.2 2000-03-30 03:07:30 edo Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
integer rtdb
|
||||
integer natom
|
||||
integer natoms
|
||||
integer nactive
|
||||
logical oactive(natom)
|
||||
logical oactive(natoms)
|
||||
integer iga_dens(2)
|
||||
double precision force(3*natom)
|
||||
double precision force(3*natoms)
|
||||
integer cntoce(*), cntobfr(*), cetobfr(*)
|
||||
double precision rdens_atom(*)
|
||||
c
|
||||
call dftgh_gridv0(rtdb, iga_dens, force, 0d0, 0, 1, natom,
|
||||
call dftgh_gridv0(rtdb, iga_dens, force, 0d0, 0, 1, natoms,
|
||||
& oactive, nactive, rdens_atom, cntoce, cntobfr,
|
||||
& cetobfr)
|
||||
c
|
||||
|
|
@ -32,19 +32,19 @@ c
|
|||
c Generalized to do XC gradient and hessian
|
||||
c
|
||||
Subroutine dftgh_gridv0(rtdb, iga_dens, force, hess, g_rhs,
|
||||
& calc_type, natom, oactive, nactive,
|
||||
& calc_type, natoms, oactive, nactive,
|
||||
& rdens_atom, cntoce, cntobfr, cetobfr)
|
||||
c
|
||||
c $Id: dftg_gridv0.F,v 1.1 2000-03-22 03:02:16 edo Exp $
|
||||
c $Id: dftg_gridv0.F,v 1.2 2000-03-30 03:07:30 edo Exp $
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
integer rtdb
|
||||
integer natom
|
||||
integer natoms
|
||||
integer nactive
|
||||
logical oactive(natom)
|
||||
logical oactive(natoms)
|
||||
integer iga_dens(2)
|
||||
double precision Force(3*natom), hess(3*natom,3*natom)
|
||||
double precision Force(3*natoms), hess(3*natoms,3*natoms)
|
||||
integer g_rhs(*) ! [input] handles for CPKS RHS mats
|
||||
integer calc_type ! [input] 1 = XC gradient
|
||||
! 2 = XC hessian
|
||||
|
|
@ -65,8 +65,8 @@ c
|
|||
c
|
||||
c
|
||||
integer cntoce(nshells_ao), cntobfr(2,nshells_ao),
|
||||
& cetobfr(2,natom)
|
||||
double precision rdens_atom(ipol*natom*natom)
|
||||
& cetobfr(2,natoms)
|
||||
double precision rdens_atom(ipol*natoms*natoms)
|
||||
logical LResult, GRAD, GRADCHI, HESSCHI, DER3CHI
|
||||
logical xc_chkgrad
|
||||
external xc_chkgrad
|
||||
|
|
@ -86,8 +86,12 @@ c
|
|||
& lPmat, iPmat, lrho, irho, lBmat, iBmat, lDmat, iDmat,
|
||||
& lhh, ihh, landex, iandex, idelrho, ldelrho,
|
||||
& ltags, itags
|
||||
logical do_force, do_hess, do_cpks_r,grid_on_file
|
||||
logical do_force, do_hess, do_cpks_r
|
||||
integer l_Amat2, k_Amat2, l_Cmat2, k_Cmat2, len_Bmat, len_Dmat
|
||||
integer lzzz, izzz, lff, iff, lppp, ippp
|
||||
integer libf_ao, iibf_ao,k_curatoms,l_curatoms
|
||||
integer lrchi_atom, irchi_atom,
|
||||
& lrdelchi_atom, irdelchi_atom
|
||||
integer nxtask
|
||||
external nxtask
|
||||
c
|
||||
|
|
@ -139,16 +143,16 @@ c write(6,*)'*** dftgh_quadv0, hesschi der3chi =',hesschi,der3chi
|
|||
c
|
||||
npol = (ipol*(ipol+1))/2
|
||||
c
|
||||
nxyz = 3*natom
|
||||
nxyz = 3*natoms
|
||||
c
|
||||
if (.not.MA_Alloc_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
|
||||
& call errquit('dftg_gridv0: cannot allocate xyz',0)
|
||||
if (.not.MA_Push_Get(MT_Dbl,natom,'charge',lcharge,icharge))
|
||||
if (.not.MA_Push_Get(MT_Dbl,natoms,'charge',lcharge,icharge))
|
||||
& call errquit('dftg_gridv0: cannot allocate charge',0)
|
||||
if (.not.MA_Push_Get(MT_Byte,natom*16,'tags',ltags,itags))
|
||||
if (.not.MA_Push_Get(MT_Byte,natoms*16,'tags',ltags,itags))
|
||||
& call errquit('dftg_gridv0: cannot allocate tags',0)
|
||||
c
|
||||
LResult = geom_cart_get(geom, natom, byte_mb(itags),
|
||||
LResult = geom_cart_get(geom, natoms, byte_mb(itags),
|
||||
& dbl_mb(ixyz), dbl_mb(icharge))
|
||||
c
|
||||
if (.not.ma_pop_stack(ltags))
|
||||
|
|
@ -159,7 +163,7 @@ c
|
|||
c get max dimens of atom block
|
||||
c
|
||||
max_at_bf = 0
|
||||
do iatom = 1, natom
|
||||
do iatom = 1, natoms
|
||||
if (.not. bas_ce2bfr(AO_bas_han, iatom, ilo, ihi))
|
||||
& call errquit('dftg_gridv0: bas_ce2bfr failed', iatom)
|
||||
max_at_bf = max(max_at_bf, ihi-ilo+1)
|
||||
|
|
@ -193,20 +197,20 @@ c
|
|||
if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',11)
|
||||
nd_qwght = 1
|
||||
if (ldew)nd_qwght = 3*n_per_rec*natom
|
||||
if (ldew)nd_qwght = 3*n_per_rec*natoms
|
||||
if (.not.MA_Push_get(mt_dbl,nd_qwght,'qwght',ld_qwght,id_qwght))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',11)
|
||||
nd_func=1
|
||||
if (ldew) nd_func=n_per_rec
|
||||
if (.not.MA_Push_get(mt_dbl,nd_func,'func',lfunc,ifunc))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',11)
|
||||
if (.not.MA_Push_get(mt_dbl,natom,'rq0',lrq0,irq0))
|
||||
if (.not.MA_Push_get(mt_dbl,natoms,'rq0',lrq0,irq0))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',12)
|
||||
if (.not.MA_Push_get(mt_dbl,natom*n_per_rec,'rq',lrq,irq))
|
||||
if (.not.MA_Push_get(mt_dbl,natoms*n_per_rec,'rq',lrq,irq))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',13)
|
||||
if (.not.MA_Push_get(mt_int,natom,'atom list',lniz,iniz))
|
||||
if (.not.MA_Push_get(mt_int,natoms,'atom list',lniz,iniz))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',13)
|
||||
if (.not.MA_Push_get(mt_int,natom,'atom nbf',lfin,ifin))
|
||||
if (.not.MA_Push_get(mt_int,natoms,'atom nbf',lfin,ifin))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',13)
|
||||
if (.not.MA_Push_Get(MT_Dbl,NXprim,'expo',lexpo,iexpo))
|
||||
& call errquit(' exiting in dftg_gridv0: insuff stack',13)
|
||||
|
|
@ -252,12 +256,34 @@ c !!! Appears that this array is never used !!! ask Jeff !!!
|
|||
c
|
||||
c stuff for neighbours list
|
||||
c
|
||||
if (.not.MA_Push_Get(MT_int,natom,'iandex',landex,iandex))
|
||||
if (.not.MA_Push_Get(MT_int,natoms,'iandex',landex,iandex))
|
||||
& call errquit(' exiting in quadvxc0: insuff stack',77)
|
||||
if (.not. rtdb_get(rtdb, 'dft:grid_on_file', mt_log, 1,
|
||||
& grid_on_file))then
|
||||
grid_on_file = .false.
|
||||
endif
|
||||
c
|
||||
c
|
||||
c allocate scratch array for grid compression
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Int, nbf_ao, 'ibf_ao',
|
||||
& libf_ao, iibf_ao))call errquit
|
||||
& (' exiting in dftg_gridv0: insuff stack',9)
|
||||
if (.not.MA_Push_Get(MT_Dbl,natoms,'rchi_atom',
|
||||
& lrchi_atom,irchi_atom))
|
||||
& call errquit('dft_scf: cannot allocate rchi_atom',0)
|
||||
if (.not.MA_Push_Get(MT_Dbl,natoms,'rdelchi_atom',
|
||||
& lrdelchi_atom,irdelchi_atom))call errquit
|
||||
& ('dftg_gridv0: cannot allocate rdelchi_atom',0)
|
||||
if (.not. MA_Push_Get(mt_dbl, nbf_ao, 'PPPP', lppp,
|
||||
& ippp))call errquit
|
||||
& (' dftg_gridv0: not enough stack',13)
|
||||
c
|
||||
c !!! BGJ - probably don't need to allocate some of this if .not.do_force
|
||||
if (.not. MA_Push_Get(mt_dbl, nbf_ao, 'F', lff, iff))
|
||||
& call errquit(' dftg_gridv0: not enough stack',13)
|
||||
if (.not. MA_Push_Get(mt_dbl, 3*nbf_ao, 'zzz', lzzz,
|
||||
& izzz))call errquit
|
||||
& (' dftg_gridv0: not enough stack',13)
|
||||
if (.not.MA_push_get(MT_INT, natoms, 'curatoms',
|
||||
& l_curatoms, k_curatoms)) call
|
||||
& errquit('dftg_gridv0a: cant allocate curatoms',1)
|
||||
c
|
||||
call dftg_gridv0a(rtdb, iga_dens,
|
||||
& ncoef_max, NXprim, max_at_bf,
|
||||
|
|
@ -271,13 +297,14 @@ c
|
|||
& dbl_mb(iBmat), dbl_mb(iDmat),
|
||||
& dbl_mb(iFmat), dbl_mb(iPmat),
|
||||
& dbl_mb(irq), dbl_mb(irq0), dbl_mb(ihh),
|
||||
, int_mb(iibf_ao),int_mb(k_curatoms),
|
||||
, dbl_mb(irchi_atom),dbl_mb(irdelchi_atom),dbl_mb(ippp),
|
||||
, dbl_mb(iff),dbl_mb(izzz),
|
||||
& force, hess, g_rhs, do_force, do_hess,
|
||||
& do_cpks_r, natom, oactive, nactive,
|
||||
& rdens_atom, cntoce, cntobfr, cetobfr)
|
||||
if (.not.ma_pop_stack(landex))
|
||||
& call errquit('dftg_gridv0: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lhh))
|
||||
& call errquit('dftg_gridv0: cannot pop stack',0)
|
||||
& do_cpks_r, natoms, oactive, nactive,
|
||||
& rdens_atom, cntobfr, cetobfr)
|
||||
if (.not.ma_chop_stack(lhh))
|
||||
& call errquit('dftg_gridv0: cannot chop stack',1)
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
if (.not.ma_pop_stack(lDmat))
|
||||
& call errquit('dftg_gridv0: cannot pop stack',0)
|
||||
|
|
@ -306,7 +333,7 @@ c
|
|||
write(LuOut,'(a,f20.10)')' integrated electron density: ',
|
||||
& rho_int
|
||||
write(LuOut,'(a,i20)')' no. of quad. points/center : ',
|
||||
& nqtot/natom
|
||||
& nqtot/natoms
|
||||
write(LuOut,*)
|
||||
end if
|
||||
c
|
||||
|
|
|
|||
|
|
@ -4,28 +4,17 @@
|
|||
& func, d_qwght, qwght, qxyz, xyz,
|
||||
& iandex, expo, ccoef, Amat, Amat2, Cmat,
|
||||
& Cmat2, Bmat, Dmat, Fmat, Pmat, rq,
|
||||
& rq0, hh, force, hess, g_rhs, do_force,
|
||||
& rq0, hh,
|
||||
. ibf_ao,curatoms,rchi_atom,rdelchi_atom,ppp,ff,zzz,
|
||||
, force, hess, g_rhs, do_force,
|
||||
& do_hess, do_cpks_r, natoms, oactive,
|
||||
& nactive, rdens_atom, cntoce, cntobfr,
|
||||
& nactive, rdens_atom, cntobfr,
|
||||
& cetobfr)
|
||||
c
|
||||
c $Id: dftg_gridv0a.F,v 1.2 2000-03-27 19:19:24 edo Exp $
|
||||
c $Id: dftg_gridv0a.F,v 1.3 2000-03-30 03:07:31 edo Exp $
|
||||
c
|
||||
implicit none
|
||||
integer rtdb
|
||||
integer natoms
|
||||
integer ncoef_max, NXprim, max_at_bf
|
||||
integer nqtot
|
||||
integer nactive
|
||||
logical oactive(natoms)
|
||||
logical grad, gradchi, hesschi, der3chi
|
||||
logical xc_chkgrad
|
||||
external xc_chkgrad
|
||||
logical ldew
|
||||
logical do_force, do_hess, do_cpks_r
|
||||
c
|
||||
#include "bas.fh"
|
||||
#include "geom.fh"
|
||||
#include "rtdb.fh"
|
||||
#include "cdft.fh"
|
||||
#include "stdio.fh"
|
||||
|
|
@ -38,9 +27,19 @@ c
|
|||
c !!! BGJ test !!!
|
||||
#include "bgj.fh"
|
||||
c !!! BGJ test !!!
|
||||
integer rtdb
|
||||
integer natoms
|
||||
integer ncoef_max, NXprim, max_at_bf
|
||||
integer nqtot
|
||||
integer nactive
|
||||
logical oactive(natoms)
|
||||
logical grad, gradchi, hesschi, der3chi
|
||||
logical xc_chkgrad
|
||||
external xc_chkgrad
|
||||
logical ldew
|
||||
logical do_force, do_hess, do_cpks_r
|
||||
c
|
||||
c
|
||||
integer cntoce(nshells_ao), cntobfr(2,nshells_ao),
|
||||
integer cntobfr(2,nshells_ao),
|
||||
& cetobfr(2,natoms)
|
||||
double precision rdens_atom(ipol*natoms*natoms)
|
||||
double precision Force(3*natoms), hess(3*natoms,3*natoms)
|
||||
|
|
@ -60,6 +59,9 @@ c
|
|||
double precision xyz(3,natoms)
|
||||
integer iandex(natoms)
|
||||
integer iniz(natoms), ifin(natoms)
|
||||
integer ibf_ao(nbf_ao),curatoms(natoms)
|
||||
double precision rchi_atom(natoms),rdelchi_atom(3,natoms)
|
||||
double precision ppp(*),ff(*),zzz(*)
|
||||
double precision ccoef(ncoef_max), expo(NXprim)
|
||||
double precision rq(n_per_rec,natoms)
|
||||
double precision rq0(natoms)
|
||||
|
|
@ -78,21 +80,11 @@ c
|
|||
c
|
||||
c local crap
|
||||
c
|
||||
integer iang, npol
|
||||
integer npol
|
||||
integer iqsh
|
||||
integer me, nproc
|
||||
integer mbf_ao, nbf_max
|
||||
integer me
|
||||
integer nbf_max
|
||||
integer NXang, llmax, MXcont, max_elem_ang_scr
|
||||
integer ibp, lbp
|
||||
integer m4
|
||||
integer nscr, iscr, lscr, icrad, iprad, lheschi_ao, iheschi_ao,
|
||||
& ld3chi_ao, id3chi_ao
|
||||
integer lzzz, izzz, lff, iff, lppp, ippp,
|
||||
& lrscr, irscr
|
||||
integer libf_ao, iibf_ao, lchi_ao, ichi_ao,
|
||||
& lgradchi_ao, idelchi_ao
|
||||
integer lrchi_atom, irchi_atom,
|
||||
& lrdelchi_atom, irdelchi_atom,icmp,lcmp
|
||||
integer nbf_ao_mxnbf_ce2
|
||||
double precision acc_AO_gauss
|
||||
double precision Ex, Ec, rho_int,rad
|
||||
|
|
@ -102,11 +94,8 @@ c
|
|||
logical grid_file_rewind
|
||||
external grid_file_rewind,grid_nbfm
|
||||
c
|
||||
integer l_curatoms, k_curatoms, ncuratoms, l_drhonuc, k_drhonuc,
|
||||
& l_ddelrhonuc, k_ddelrhonuc, l_d2scr, k_d2scr
|
||||
c
|
||||
me = ga_nodeid()
|
||||
nproc = ga_nnodes()
|
||||
Ec = 0.d0
|
||||
Ex = 0.d0
|
||||
if (bgj_print() .gt. 0)
|
||||
|
|
@ -120,11 +109,18 @@ c Define the number of basis functions in the AO Gaussian basis set.
|
|||
c
|
||||
if (.not. bas_ncontr_cn_max(ao_bas_han, MXcont))
|
||||
& call errquit('Exiting in dftg_quadv0.',1)
|
||||
c
|
||||
|
||||
if (.not. bas_high_angular(AO_bas_han, NXang))
|
||||
& call errquit('exiting in dftg_quadv0.',7)
|
||||
llMAX=((NXang+2)*(NXang+1))/2
|
||||
c scratch for eval_func
|
||||
if (NXang .ge. 3)then
|
||||
max_elem_ang_scr = 28
|
||||
elseif (NXang .ge. 1)then
|
||||
max_elem_ang_scr = 9
|
||||
else
|
||||
max_elem_ang_scr = 3
|
||||
endif
|
||||
max_elem_ang_scr=max_elem_ang_scr*llMAX
|
||||
c
|
||||
if (.not.bas_nbf_cn_max(AO_bas_han, nbf_max))
|
||||
& call errquit('exiting in dftg_quadv0.',7)
|
||||
|
|
@ -134,7 +130,6 @@ c Define the max number of basis functions, contractions,
|
|||
c and primitives comprising any given shell of this basis set.
|
||||
c
|
||||
acc_AO_gauss= iAOacc
|
||||
c
|
||||
c
|
||||
if (.not. rtdb_get(rtdb, 'dft:weight derivatives', mt_log, 1,
|
||||
& ldew))ldew=.true.
|
||||
|
|
@ -148,24 +143,12 @@ c
|
|||
& ' derivative weights not available for Delley weights '
|
||||
ldew=.false.
|
||||
endif
|
||||
if (.not. rtdb_get(rtdb, 'dft:ssw', mt_log, 1, lssw))
|
||||
& lssw=.false.
|
||||
if(lssw) then
|
||||
if (.not.rtdb_cget(rtdb, 'dft:ssweights', 1, whichssw))
|
||||
& call errquit('dftgq0: rtdb_cget failed', 1700)
|
||||
endif
|
||||
c
|
||||
c see if any Gradient Corrected Functional is not null
|
||||
c
|
||||
grad=xc_chkgrad(xfac,cfac)
|
||||
c
|
||||
c allocate scratch array for grid compression
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Int, nbf_ao, 'ibf_ao',
|
||||
& libf_ao, iibf_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
if (.not.MA_Push_Get(mt_int,n_per_rec,'compr scr',lcmp,icmp))
|
||||
& call errquit('dftg_quadv0_a: cannot allocate compr',0)
|
||||
c
|
||||
if (grad)then
|
||||
gradchi = .true.
|
||||
|
|
@ -179,35 +162,6 @@ c
|
|||
c write(6,*)'*** dftg_quadv0a: hesschi der3chi',hesschi,der3chi
|
||||
c
|
||||
npol = (ipol*(ipol+1))/2
|
||||
c scratch for eval_func
|
||||
if (NXang .ge. 3)then
|
||||
max_elem_ang_scr = 28
|
||||
elseif (NXang .ge. 1)then
|
||||
max_elem_ang_scr = 9
|
||||
else
|
||||
max_elem_ang_scr = 3
|
||||
endif
|
||||
max_elem_ang_scr=max_elem_ang_scr*llMAX
|
||||
c
|
||||
c Allocate for reduced basis function quantities early since they
|
||||
c are used in the 2nd derivative code further down
|
||||
c
|
||||
if (.not.MA_Push_Get(MT_Dbl,natoms,'rchi_atom',
|
||||
& lrchi_atom,irchi_atom))
|
||||
& call errquit('dft_scf: cannot allocate rchi_atom',0)
|
||||
if (.not.MA_Push_Get(MT_Dbl,natoms,'rdelchi_atom',
|
||||
& lrdelchi_atom,irdelchi_atom))call errquit
|
||||
& ('dft_scf: cannot allocate rdelchi_atom',0)
|
||||
if (.not. MA_Push_Get(mt_dbl, nbf_ao, 'PPPP', lppp,
|
||||
& ippp))call errquit
|
||||
& (' dftg_quadv0: not enough stack',13)
|
||||
c
|
||||
c !!! BGJ - probably don't need to allocate some of this if .not.do_force
|
||||
if (.not. MA_Push_Get(mt_dbl, nbf_ao, 'F', lff, iff))
|
||||
& call errquit(' dftg_quadv0: not enough stack',13)
|
||||
if (.not. MA_Push_Get(mt_dbl, 3*nbf_ao, 'zzz', lzzz,
|
||||
& izzz))call errquit
|
||||
& (' dftg_quadv0: not enough stack',13)
|
||||
c
|
||||
c rewind grid pts file
|
||||
c
|
||||
|
|
@ -236,379 +190,26 @@ c
|
|||
call grid_repack(dbl_mb(irqbuf), qxyz, qwght, nqpts,
|
||||
, rad,istep)
|
||||
if(nqpts.eq.0) goto 110
|
||||
|
||||
c Evaluate the AO basis set at each of the quad. points.
|
||||
c allocate arrays for exponents and contraction coefficients
|
||||
c within int_eval_gbsets mbf is # of orbs with nonzero
|
||||
c points in sampling set. iniz is a mapping array.
|
||||
c currently set up only for cartesian (6d) sets
|
||||
c
|
||||
mbf_ao = nbf_ao
|
||||
call ifill(mbf_ao, 0, int_mb(iibf_ao), 1)
|
||||
mbf_ao=grid_nbfm( AO_bas_han, natoms,
|
||||
& ictr_buf,rad,xyz,qxyz,nqpts,
|
||||
& int_mb(iibf_ao), iniz, ifin, expo,
|
||||
. minexp,ldiff, acc_AO_gauss)
|
||||
if (mbf_ao .eq. 0)goto 110
|
||||
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
c
|
||||
c Set up space for derivatives of density wrt significant nuclei
|
||||
c
|
||||
c write(*,*)'allocating deriv density space'
|
||||
if (.not.MA_push_get(MT_INT, natoms, 'curatoms',
|
||||
& l_curatoms, k_curatoms)) call
|
||||
& errquit('dftg_quadv0a: cant allocate curatoms',1)
|
||||
if (.not.MA_push_get(MT_DBL, ipol*nqpts*3*nactive,
|
||||
& 'drhonuc', l_drhonuc, k_drhonuc)) call
|
||||
& errquit('dftg_quadv0a: cannot allocate drhonuc',1)
|
||||
if (grad) then
|
||||
if (.not.MA_push_get(MT_DBL,
|
||||
& ipol*3*nqpts*3*nactive, 'ddelrhonuc',
|
||||
& l_ddelrhonuc, k_ddelrhonuc)) call
|
||||
& errquit('dftg_quadv0a: cant alloc ddelrhonuc',
|
||||
& 1)
|
||||
endif
|
||||
endif
|
||||
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Dbl, nqpts*mbf_ao, 'chi_ao',
|
||||
& lchi_ao, ichi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
if (.not. MA_Push_Get(MT_Dbl, 3*nqpts*mbf_ao,
|
||||
& 'gradchi_ao', lgradchi_ao,idelchi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
c
|
||||
c hessian is used only for derivatives; set array to length 1
|
||||
c here
|
||||
c
|
||||
c
|
||||
if (hesschi)then
|
||||
m4 = 6*nqpts*mbf_ao
|
||||
else
|
||||
m4 = 1
|
||||
endif
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Dbl, m4, 'heschi_ao',
|
||||
& lheschi_ao, iheschi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
c
|
||||
if (der3chi)then
|
||||
m4 = 10*nqpts*mbf_ao
|
||||
else
|
||||
m4 = 1
|
||||
endif
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Dbl, m4, 'd3chi_ao',
|
||||
& ld3chi_ao, id3chi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',19)
|
||||
c
|
||||
c compute scratch space for int_eval_gbsets
|
||||
c
|
||||
if (der3chi) then
|
||||
c
|
||||
c Allocate scratch space.
|
||||
c mxcont is used here instead of NXprim as in new_eval_gbsets
|
||||
c this is correct, but sometimes fails, hence the choice
|
||||
c in new_eval_gsets. here, the length is also doubled...
|
||||
c
|
||||
icrad = 1 ! contracted radial parts
|
||||
iprad = icrad + 4*nqpts*MXcont ! primitive radial parts
|
||||
iang = iprad ! also angular parts
|
||||
c
|
||||
c note 4*NXprim for longer primitive and contracted radial
|
||||
c scratch storage
|
||||
c
|
||||
nscr=iang+max(4*nqpts*NXprim,max_elem_ang_scr*nqpts)
|
||||
c
|
||||
else if (hesschi)then
|
||||
c
|
||||
c Allocate scratch space.
|
||||
c mxcont is used here instead of NXprim as in new_eval_gbsets
|
||||
c this is correct, but sometimes fails, hence the choice
|
||||
c in new_eval_gsets. here, the length is also doubled...
|
||||
c
|
||||
icrad = 1 ! contracted radial parts
|
||||
iprad = icrad + 3*nqpts*MXcont ! primitive radial parts
|
||||
iang = iprad ! also angular parts
|
||||
c
|
||||
c note 3*NXprim for longer primitive and contracted radial
|
||||
c scratch storage
|
||||
c
|
||||
nscr=iang+max(3*nqpts*NXprim,max_elem_ang_scr*nqpts)
|
||||
c
|
||||
elseif (gradchi)then
|
||||
c
|
||||
c Allocate scratch space.
|
||||
c
|
||||
icrad = 1
|
||||
iprad = icrad + 2*nqpts*NXprim
|
||||
iang = iprad
|
||||
nscr = iang+max(2*nqpts*NXprim,9*nqpts*llMAX) - 1
|
||||
nscr = 2*nscr
|
||||
else
|
||||
c
|
||||
c Allocate scratch space.
|
||||
|
||||
icrad = 1
|
||||
iprad = icrad + nqpts*NXprim
|
||||
iang = iprad
|
||||
nscr = iang + max(nqpts*NXprim,6*nqpts*llMAX)
|
||||
& - 1
|
||||
nscr = 2*nscr
|
||||
c
|
||||
endif
|
||||
c
|
||||
if (.not.MA_Push_Get(mt_dbl, nscr, 'scrCD', lscr, iscr))
|
||||
& call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',12)
|
||||
c
|
||||
call qdist(rq0, rq, qxyz, xyz, nqpts, natoms)
|
||||
if (der3chi) then
|
||||
call new_eval_gbsets_d3(ao_bas_han,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), dbl_mb(id3chi_ao),
|
||||
& int_mb(iibf_ao), 0, rq0, rq,
|
||||
& qxyz, xyz, nqpts, nqpts, mbf_ao,
|
||||
& natoms, dbl_mb(iscr), nscr, expo, NXprim,
|
||||
& ccoef, ncoef_max,
|
||||
& iniz,ifin)
|
||||
else
|
||||
call new_eval_gbsets(ao_bas_han, gradchi,
|
||||
& hesschi, dbl_mb(ichi_ao),
|
||||
& dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao),
|
||||
& int_mb(iibf_ao), 0, rq0,
|
||||
& rq, qxyz, xyz, nqpts,
|
||||
& nqpts, mbf_ao, natoms,
|
||||
. dbl_mb(iscr),
|
||||
& nscr, expo, NXprim, ccoef,
|
||||
& ncoef_max,
|
||||
& iniz,ifin)
|
||||
endif
|
||||
if (.not.ma_pop_stack(lscr))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
c
|
||||
c Pre-compute reduced Xi(r) and dXi(r) over atoms
|
||||
c reduce chi(nq,mbf_ao) to chi(natoms)
|
||||
c reduce delchi(nq,3,mbf_ao) to delchi(natoms)
|
||||
call util_rmsatbf(nqpts, natoms,iniz,ifin,
|
||||
, dbl_mb(ichi_ao),dbl_mb(irchi_atom))
|
||||
if(grad)
|
||||
. call util_rmsatbf(nqpts, natoms,iniz,ifin,
|
||||
, dbl_mb(idelchi_ao),dbl_mb(irdelchi_atom))
|
||||
! call grid_rchi3(ntot_cube+1,
|
||||
! , nqpts, natoms,iniz,ifin,
|
||||
! , dbl_mb(ichi_ao),dbl_mb(idelchi_ao),
|
||||
! , dbl_mb(irchi_atom),dbl_mb(irdelchi_atom),
|
||||
! , grad)
|
||||
c
|
||||
if (.not.MA_Push_Get(mt_dbl, nqpts*natoms, 'bp',
|
||||
& lbp, ibp))call errquit
|
||||
& (' exiting in quadvxc0: insuff stack',7)
|
||||
|
||||
if (.not.MA_Push_Get(mt_dbl,nqpts,'rho scr',
|
||||
& lrscr,irscr))call errquit
|
||||
& ('dftg_quadv0a: cannot allocate rho scr',12)
|
||||
c
|
||||
call xc_rhogen(tol_rho, AO_bas_han, iga_dens, max_at_bf,
|
||||
& natoms, ipol, nqpts, nbf_ao, mbf_ao,
|
||||
& ldelley, grad, Fmat, Pmat, Bmat,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& int_mb(iibf_ao), iniz, ifin, rho, delrho,
|
||||
& dbl_mb(ibp), dbl_mb(irscr),
|
||||
& dbl_mb(irchi_atom), dbl_mb(irdelchi_atom),
|
||||
& rdens_atom, cetobfr)
|
||||
c
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
call oactive_to_curatoms(oactive, natoms,
|
||||
& int_mb(k_curatoms), ncuratoms)
|
||||
call dfill(ipol*nqpts*3*nactive,0.D0,
|
||||
& dbl_mb(k_drhonuc),1)
|
||||
if (grad)
|
||||
& call dfill(ipol*3*nqpts*3*nactive,0.D0,
|
||||
& dbl_mb(k_ddelrhonuc),1)
|
||||
c write(*,*)'calling xc_rhogen_nucder'
|
||||
#if 1
|
||||
call xc_rhogen_nucder(tol_rho, AO_bas_han, iga_dens,
|
||||
& max_at_bf, natoms, int_mb(k_curatoms), ncuratoms,
|
||||
& ipol, nqpts, nbf_ao, mbf_ao, grad,
|
||||
& Fmat, Pmat, Bmat, Dmat,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), int_mb(iibf_ao), iniz, ifin,
|
||||
& dbl_mb(k_drhonuc), dbl_mb(k_ddelrhonuc),
|
||||
& dbl_mb(irscr), dbl_mb(irchi_atom),
|
||||
& dbl_mb(irdelchi_atom), rdens_atom, cetobfr)
|
||||
#endif
|
||||
endif
|
||||
c
|
||||
if (.not.ma_pop_stack(lrscr))
|
||||
& call errquit('xc_quadv0_a: cannot pop stack',0)
|
||||
c
|
||||
if (.not.ma_pop_stack(lbp))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
c Evaluate the integrated charge density.
|
||||
c
|
||||
nqtot = nqtot + nqpts
|
||||
call dftg_gridv0b(nqpts,rad,ictr_buf,iga_dens,
|
||||
& ncoef_max, NXprim, max_at_bf,
|
||||
& rho_int, nqtot, iniz, ifin, rho, delrho,
|
||||
& func, d_qwght, qwght, qxyz, xyz,
|
||||
& iandex, expo, ccoef, Amat, Amat2, Cmat,
|
||||
& Cmat2, Bmat, Dmat, Fmat, Pmat, rq,
|
||||
& rq0, hh,
|
||||
. ibf_ao,curatoms,rchi_atom,rdelchi_atom,ppp,ff,zzz,
|
||||
, force, hess, g_rhs, do_force,
|
||||
, llmax,max_elem_ang_scr,nbf_ao_mxnbf_ce2,mxcont,
|
||||
, acc_ao_gauss, ldew,grad,gradchi,hesschi,der3chi,npol,
|
||||
& do_hess, do_cpks_r, natoms, oactive,
|
||||
& nactive, rdens_atom,
|
||||
& cetobfr)
|
||||
|
||||
rho_int = rho_int +
|
||||
& ddot(nqpts, qwght, 1, rho, 1)
|
||||
|
||||
call xc_eval_fnl(rho, delrho, Amat, Amat2, Cmat, Cmat2,
|
||||
& nqpts, Ex, Ec, qwght, GRAD, ldew, func,
|
||||
& do_hess .or. do_cpks_r)
|
||||
c
|
||||
c
|
||||
c Exchange-correlation contribution
|
||||
c
|
||||
c exchange-correlation contribution, with or without gradient
|
||||
c terms
|
||||
c
|
||||
if (do_force) then
|
||||
c
|
||||
call dftg_quadv0_b(ictr_buf, ldew, AO_bas_han, ipol,
|
||||
& grad, nbf_ao, mbf_ao, natoms, nqpts,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), Amat, Cmat,
|
||||
& int_mb(iibf_ao), iniz, ifin,
|
||||
& dbl_mb(ippp), dbl_mb(iff),
|
||||
& dbl_mb(izzz), hh, iga_dens, force,
|
||||
& oactive, nactive)
|
||||
c
|
||||
else if (do_hess) then
|
||||
c
|
||||
if (.not.MA_push_get(mt_dbl, 15*nqpts, 'd2scr',
|
||||
& l_d2scr, k_d2scr)) call errquit
|
||||
& ('dftg_quadv0: not enough stack for d2scr',0)
|
||||
c
|
||||
c Note that k_d2scr is passed in three times, since it is also
|
||||
c used for the scratch arrays ff and s. This is ok since by the
|
||||
c time ff and s are used, the routine is done with the array scr.
|
||||
c -- what about using Bmat? !!!
|
||||
call xc_d2expl(tol_rho, dbl_mb(k_d2scr), Amat, Amat2,
|
||||
& Bmat, Cmat, Cmat2, Dmat, Fmat, Pmat,
|
||||
& dbl_mb(k_d2scr), dbl_mb(k_d2scr+3*nqpts),
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), dbl_mb(id3chi_ao),
|
||||
& int_mb(k_curatoms), ncuratoms, ipol, nqpts,
|
||||
& nbf_ao, mbf_ao, grad, AO_bas_han, natoms, iniz,
|
||||
& ifin, dbl_mb(k_drhonuc), dbl_mb(k_ddelrhonuc),
|
||||
& delrho, iga_dens, hess, int_mb(iibf_ao),
|
||||
& dbl_mb(irchi_atom), dbl_mb(irdelchi_atom),
|
||||
& rdens_atom, cetobfr)
|
||||
if (.not.ma_pop_stack(l_d2scr))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
else if (do_cpks_r) then
|
||||
c
|
||||
c The CPKS RHS is conveniently done in two parts: the part
|
||||
c involving the derivs of the density looks just like the CPKS
|
||||
c LHS, so that code is called for that part, while the part
|
||||
c involving only the derivs of the basis functions is calculated
|
||||
c with a separate routine
|
||||
c
|
||||
c First part: LHS-like contribution
|
||||
c
|
||||
c Combine derivs of functional with deriv densities to get coeffs
|
||||
c of bfn pairs in CPKS LHS - these are written back into rho and
|
||||
c delrho
|
||||
c
|
||||
call xc_cpks_coeff(dbl_mb(k_drhonuc),
|
||||
& dbl_mb(k_ddelrhonuc), Amat2, Cmat2, Cmat, delrho,
|
||||
& 3*ncuratoms, ipol, nqpts, GRAD)
|
||||
#if 1
|
||||
c write(*,*)'cpks first part'
|
||||
call xc_tabcd_cpks_lhs(tol_rho, Fmat, Pmat,
|
||||
& dbl_mb(k_drhonuc), Bmat, dbl_mb(k_ddelrhonuc),
|
||||
& Dmat, Dbl_MB(ichi_ao), Dbl_MB(idelchi_ao),
|
||||
& 3*ncuratoms, ipol, nqpts, mbf_ao,
|
||||
& nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2, GRAD,
|
||||
& AO_bas_han, natoms, iniz, ifin, g_rhs,
|
||||
& INT_MB(iibf_ao), dbl_mb(irchi_atom),
|
||||
& dbl_mb(irdelchi_atom), rdens_atom, cetobfr)
|
||||
c write(*,*)'Cmat after lhs'
|
||||
c call output(cmat,1,nq,1,6,nq,6,1)
|
||||
#endif
|
||||
#if 1
|
||||
c
|
||||
c Second part: Basis function derivative contribution
|
||||
c
|
||||
c write(*,*)'cpks second part'
|
||||
c
|
||||
c Transform Cmat to delrho form for use in xc_cpks_tabcd_rhs
|
||||
c
|
||||
if (grad) then
|
||||
call transform_Cmat(delrho, Cmat, ipol, nqpts)
|
||||
endif
|
||||
call xc_tabcd_cpks_rhs(tol_rho, Fmat, Pmat, Amat,
|
||||
& Bmat, Cmat, Dmat, Dbl_MB(ichi_ao),
|
||||
& Dbl_MB(idelchi_ao), dbl_mb(iheschi_ao),
|
||||
& int_mb(k_curatoms),
|
||||
& ncuratoms, ipol, nqpts, mbf_ao,
|
||||
& nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2, GRAD,
|
||||
& AO_bas_han, natoms, iniz, ifin, g_rhs,
|
||||
& INT_MB(iibf_ao), dbl_mb(irchi_atom),
|
||||
& dbl_mb(irdelchi_atom), rdens_atom, cetobfr)
|
||||
#endif
|
||||
c
|
||||
endif
|
||||
c
|
||||
if (.not.ma_pop_stack(ld3chi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lheschi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lgradchi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lchi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
if (grad) then
|
||||
if (.not.ma_pop_stack(l_ddelrhonuc))
|
||||
& call errquit('xc_quadv0_a: cannot pop ddelrhonuc',
|
||||
& 0)
|
||||
endif
|
||||
if (.not.ma_pop_stack(l_drhonuc))
|
||||
& call errquit('xc_quadv0_a: cannot pop drhonuc',0)
|
||||
if (.not.ma_pop_stack(l_curatoms))
|
||||
& call errquit('xc_quadv0_a: cannot pop curatoms',0)
|
||||
endif
|
||||
c
|
||||
c
|
||||
cedo no weight derivatives for the moment
|
||||
cedo if (do_force.and.ldew)then
|
||||
cedo write(0,*) ictr_buf,natoms,nqpts
|
||||
cedo do i=1,natoms
|
||||
cedo iandex(ii)=ii
|
||||
cedo enddo
|
||||
cedo call dftg_wderiv(ictr_buf, natoms, nqpts, nqpts,
|
||||
cedo . iandex,
|
||||
cedo & qcenters, func, d_qwght,
|
||||
cedo & force, oactive, nactive)
|
||||
cedo endif
|
||||
c
|
||||
c
|
||||
110 ntot_cube=ntot_cube+1
|
||||
200 continue ! bottom of shell loop
|
||||
if (.not.ma_pop_stack(lrqbuf))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack', 0)
|
||||
if (.not.ma_pop_stack(lzzz))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lff))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lppp))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lrdelchi_atom))
|
||||
& call errquit('xc_quadv0_a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lrchi_atom))
|
||||
& call errquit('xc_quadv0_a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lcmp))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack', 0)
|
||||
c
|
||||
if (.not.ma_pop_stack(libf_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
return
|
||||
c
|
||||
end
|
||||
|
|
|
|||
455
src/dftgrad/dftg_gridv0b.F
Normal file
455
src/dftgrad/dftg_gridv0b.F
Normal file
|
|
@ -0,0 +1,455 @@
|
|||
Subroutine dftg_gridv0b(nqpts,rad,ictr_buf,iga_dens,
|
||||
& ncoef_max, NXprim, max_at_bf,
|
||||
& rho_int, nqtot, iniz, ifin, rho, delrho,
|
||||
& func, d_qwght, qwght, qxyz, xyz,
|
||||
& iandex, expo, ccoef, Amat, Amat2, Cmat,
|
||||
& Cmat2, Bmat, Dmat, Fmat, Pmat, rq,
|
||||
& rq0, hh,
|
||||
. ibf_ao,curatoms,rchi_atom,rdelchi_atom,ppp,ff,zzz,
|
||||
, force, hess, g_rhs, do_force,
|
||||
, llmax,max_elem_ang_scr,nbf_ao_mxnbf_ce2,mxcont,
|
||||
, acc_ao_gauss, ldew,grad,gradchi,hesschi,der3chi,npol,
|
||||
& do_hess, do_cpks_r, natoms, oactive,
|
||||
& nactive, rdens_atom,
|
||||
& cetobfr)
|
||||
c
|
||||
c $Id: dftg_gridv0b.F,v 1.1 2000-03-30 03:07:31 edo Exp $
|
||||
c
|
||||
implicit none
|
||||
#include "bas.fh"
|
||||
#include "cdft.fh"
|
||||
#include "stdio.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "global.fh"
|
||||
#include "msgids.fh"
|
||||
#include "cgridfile.fh"
|
||||
c
|
||||
#include "dft2drv.fh"
|
||||
c !!! BGJ test !!!
|
||||
#include "bgj.fh"
|
||||
c !!! BGJ test !!!
|
||||
integer nqpts, ictr_buf ! [in]
|
||||
double precision rad ! [in]
|
||||
integer natoms
|
||||
integer ncoef_max, NXprim, max_at_bf
|
||||
integer nqtot
|
||||
integer nactive
|
||||
logical oactive(natoms)
|
||||
logical grad, gradchi, hesschi, der3chi
|
||||
logical xc_chkgrad
|
||||
external xc_chkgrad
|
||||
logical ldew
|
||||
logical do_force, do_hess, do_cpks_r
|
||||
c
|
||||
integer cetobfr(2,natoms)
|
||||
double precision rdens_atom(ipol*natoms*natoms)
|
||||
double precision Force(3*natoms), hess(3*natoms,3*natoms)
|
||||
integer g_rhs(*)
|
||||
double precision hh(3,n_per_rec)
|
||||
double precision rho(n_per_rec,ipol*(ipol+1)/2)
|
||||
double precision delrho(n_per_rec,3,ipol)
|
||||
double precision Amat(n_per_rec,ipol), Cmat(n_per_rec,3,ipol),
|
||||
& Bmat(n_per_rec,*), Dmat(n_per_rec,3,max_at_bf)
|
||||
double precision Amat2(n_per_rec,NCOL_AMAT2)
|
||||
double precision Cmat2(n_per_rec,NCOL_CMAT2)
|
||||
double precision Pmat(max_at_bf*max_at_bf),
|
||||
& Fmat(max_at_bf*max_at_bf)
|
||||
double precision func(n_per_rec)
|
||||
double precision d_qwght(3,n_per_rec,*)
|
||||
double precision qwght(nqpts), qxyz(3,nqpts)
|
||||
double precision xyz(3,natoms)
|
||||
integer iandex(natoms)
|
||||
integer iniz(natoms), ifin(natoms)
|
||||
integer ibf_ao(nbf_ao),curatoms(natoms)
|
||||
double precision rchi_atom(natoms),rdelchi_atom(3,natoms)
|
||||
double precision ppp(*),ff(*),zzz(*)
|
||||
double precision ccoef(ncoef_max), expo(NXprim)
|
||||
double precision rq(n_per_rec,natoms)
|
||||
double precision rq0(natoms)
|
||||
c
|
||||
c external
|
||||
c
|
||||
double precision ddot
|
||||
integer nxtask
|
||||
External ddot, nxtask
|
||||
double precision eps
|
||||
parameter(eps=1.d-8)
|
||||
c
|
||||
c---> Matrix Elements of the XC Potential & Energy
|
||||
c
|
||||
integer iga_dens(2)
|
||||
c
|
||||
c local crap
|
||||
c
|
||||
integer iang, npol
|
||||
integer mbf_ao
|
||||
integer llmax, MXcont, max_elem_ang_scr
|
||||
integer ibp, lbp
|
||||
integer m4
|
||||
integer nscr, iscr, lscr, icrad, iprad, lheschi_ao, iheschi_ao,
|
||||
& ld3chi_ao, id3chi_ao
|
||||
integer lrscr, irscr
|
||||
integer lchi_ao, ichi_ao,
|
||||
& lgradchi_ao, idelchi_ao
|
||||
integer nbf_ao_mxnbf_ce2
|
||||
double precision acc_AO_gauss
|
||||
double precision Ex, Ec
|
||||
double precision rho_int ! [in/out]
|
||||
integer grid_nbfm
|
||||
external grid_nbfm
|
||||
c
|
||||
integer ncuratoms, l_drhonuc, k_drhonuc,
|
||||
& l_ddelrhonuc, k_ddelrhonuc, l_d2scr, k_d2scr
|
||||
c
|
||||
|
||||
c Evaluate the AO basis set at each of the quad. points.
|
||||
c allocate arrays for exponents and contraction coefficients
|
||||
c within int_eval_gbsets mbf is # of orbs with nonzero
|
||||
c points in sampling set. iniz is a mapping array.
|
||||
c currently set up only for cartesian (6d) sets
|
||||
c
|
||||
mbf_ao = nbf_ao
|
||||
call ifill(mbf_ao, 0, ibf_ao, 1)
|
||||
mbf_ao=grid_nbfm( AO_bas_han, natoms,
|
||||
& ictr_buf,rad,xyz,qxyz,nqpts,
|
||||
& ibf_ao, iniz, ifin, expo,
|
||||
. minexp,ldiff, acc_AO_gauss)
|
||||
if (mbf_ao .eq. 0) return
|
||||
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
c
|
||||
c Set up space for derivatives of density wrt significant nuclei
|
||||
c
|
||||
c write(*,*)'allocating deriv density space'
|
||||
if (.not.MA_push_get(MT_DBL, ipol*nqpts*3*nactive,
|
||||
& 'drhonuc', l_drhonuc, k_drhonuc)) call
|
||||
& errquit('dftg_quadv0a: cannot allocate drhonuc',1)
|
||||
if (grad) then
|
||||
if (.not.MA_push_get(MT_DBL,
|
||||
& ipol*3*nqpts*3*nactive, 'ddelrhonuc',
|
||||
& l_ddelrhonuc, k_ddelrhonuc)) call
|
||||
& errquit('dftg_quadv0a: cant alloc ddelrhonuc',
|
||||
& 1)
|
||||
endif
|
||||
endif
|
||||
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Dbl, nqpts*mbf_ao, 'chi_ao',
|
||||
& lchi_ao, ichi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
if (.not. MA_Push_Get(MT_Dbl, 3*nqpts*mbf_ao,
|
||||
& 'gradchi_ao', lgradchi_ao,idelchi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
c
|
||||
c hessian is used only for derivatives; set array to length 1
|
||||
c here
|
||||
c
|
||||
c
|
||||
if (hesschi)then
|
||||
m4 = 6*nqpts*mbf_ao
|
||||
else
|
||||
m4 = 1
|
||||
endif
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Dbl, m4, 'heschi_ao',
|
||||
& lheschi_ao, iheschi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',9)
|
||||
c
|
||||
if (der3chi)then
|
||||
m4 = 10*nqpts*mbf_ao
|
||||
else
|
||||
m4 = 1
|
||||
endif
|
||||
c
|
||||
if (.not. MA_Push_Get(MT_Dbl, m4, 'd3chi_ao',
|
||||
& ld3chi_ao, id3chi_ao))call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',19)
|
||||
c
|
||||
c compute scratch space for int_eval_gbsets
|
||||
c
|
||||
if (der3chi) then
|
||||
c
|
||||
c Allocate scratch space.
|
||||
c mxcont is used here instead of NXprim as in new_eval_gbsets
|
||||
c this is correct, but sometimes fails, hence the choice
|
||||
c in new_eval_gsets. here, the length is also doubled...
|
||||
c
|
||||
icrad = 1 ! contracted radial parts
|
||||
iprad = icrad + 4*nqpts*MXcont ! primitive radial parts
|
||||
iang = iprad ! also angular parts
|
||||
c
|
||||
c note 4*NXprim for longer primitive and contracted radial
|
||||
c scratch storage
|
||||
c
|
||||
nscr=iang+max(4*nqpts*NXprim,max_elem_ang_scr*nqpts)
|
||||
c
|
||||
else if (hesschi)then
|
||||
c
|
||||
c Allocate scratch space.
|
||||
c mxcont is used here instead of NXprim as in new_eval_gbsets
|
||||
c this is correct, but sometimes fails, hence the choice
|
||||
c in new_eval_gsets. here, the length is also doubled...
|
||||
c
|
||||
icrad = 1 ! contracted radial parts
|
||||
iprad = icrad + 3*nqpts*MXcont ! primitive radial parts
|
||||
iang = iprad ! also angular parts
|
||||
c
|
||||
c note 3*NXprim for longer primitive and contracted radial
|
||||
c scratch storage
|
||||
c
|
||||
nscr=iang+max(3*nqpts*NXprim,max_elem_ang_scr*nqpts)
|
||||
c
|
||||
elseif (gradchi)then
|
||||
c
|
||||
c Allocate scratch space.
|
||||
c
|
||||
icrad = 1
|
||||
iprad = icrad + 2*nqpts*NXprim
|
||||
iang = iprad
|
||||
nscr = iang+max(2*nqpts*NXprim,9*nqpts*llMAX) - 1
|
||||
nscr = 2*nscr
|
||||
else
|
||||
c
|
||||
c Allocate scratch space.
|
||||
|
||||
icrad = 1
|
||||
iprad = icrad + nqpts*NXprim
|
||||
iang = iprad
|
||||
nscr = iang + max(nqpts*NXprim,6*nqpts*llMAX)
|
||||
& - 1
|
||||
nscr = 2*nscr
|
||||
c
|
||||
endif
|
||||
c
|
||||
if (.not.MA_Push_Get(mt_dbl, nscr, 'scrCD', lscr, iscr))
|
||||
& call errquit
|
||||
& (' exiting in dftg_quadv0: insuff stack',12)
|
||||
c
|
||||
call qdist(rq0, rq, qxyz, xyz, nqpts, natoms)
|
||||
if (der3chi) then
|
||||
call new_eval_gbsets_d3(ao_bas_han,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), dbl_mb(id3chi_ao),
|
||||
& ibf_ao, 0, rq0, rq,
|
||||
& qxyz, xyz, nqpts, nqpts, mbf_ao,
|
||||
& natoms, dbl_mb(iscr), nscr, expo, NXprim,
|
||||
& ccoef, ncoef_max,
|
||||
& iniz,ifin)
|
||||
else
|
||||
call new_eval_gbsets(ao_bas_han, gradchi,
|
||||
& hesschi, dbl_mb(ichi_ao),
|
||||
& dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao),
|
||||
& ibf_ao, 0, rq0,
|
||||
& rq, qxyz, xyz, nqpts,
|
||||
& nqpts, mbf_ao, natoms,
|
||||
. dbl_mb(iscr),
|
||||
& nscr, expo, NXprim, ccoef,
|
||||
& ncoef_max,
|
||||
& iniz,ifin)
|
||||
endif
|
||||
if (.not.ma_pop_stack(lscr))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
c
|
||||
c Pre-compute reduced Xi(r) and dXi(r) over atoms
|
||||
c reduce chi(nq,mbf_ao) to chi(natoms)
|
||||
c reduce delchi(nq,3,mbf_ao) to delchi(natoms)
|
||||
call util_rmsatbf(nqpts, natoms,iniz,ifin,
|
||||
, dbl_mb(ichi_ao),rchi_atom)
|
||||
if(grad)
|
||||
. call util_rmsatbf(nqpts, natoms,iniz,ifin,
|
||||
, dbl_mb(idelchi_ao),rdelchi_atom)
|
||||
! call grid_rchi3(ntot_cube+1,
|
||||
! , nqpts, natoms,iniz,ifin,
|
||||
! , dbl_mb(ichi_ao),dbl_mb(idelchi_ao),
|
||||
! , rchi_atom,dbl_mb(irdelchi_atom),
|
||||
! , grad)
|
||||
c
|
||||
if (.not.MA_Push_Get(mt_dbl, nqpts*natoms, 'bp',
|
||||
& lbp, ibp))call errquit
|
||||
& (' exiting in quadvxc0: insuff stack',7)
|
||||
|
||||
if (.not.MA_Push_Get(mt_dbl,nqpts,'rho scr',
|
||||
& lrscr,irscr))call errquit
|
||||
& ('dftg_quadv0a: cannot allocate rho scr',12)
|
||||
c
|
||||
call xc_rhogen(tol_rho, AO_bas_han, iga_dens, max_at_bf,
|
||||
& natoms, ipol, nqpts, nbf_ao, mbf_ao,
|
||||
& ldelley, grad, Fmat, Pmat, Bmat,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& ibf_ao, iniz, ifin, rho, delrho,
|
||||
& dbl_mb(ibp), dbl_mb(irscr),
|
||||
& rchi_atom, rdelchi_atom,
|
||||
& rdens_atom, cetobfr)
|
||||
c
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
call oactive_to_curatoms(oactive, natoms,
|
||||
& curatoms, ncuratoms)
|
||||
call dfill(ipol*nqpts*3*nactive,0.D0,
|
||||
& dbl_mb(k_drhonuc),1)
|
||||
if (grad)
|
||||
& call dfill(ipol*3*nqpts*3*nactive,0.D0,
|
||||
& dbl_mb(k_ddelrhonuc),1)
|
||||
c write(*,*)'calling xc_rhogen_nucder'
|
||||
#if 1
|
||||
call xc_rhogen_nucder(tol_rho, AO_bas_han, iga_dens,
|
||||
& max_at_bf, natoms, curatoms, ncuratoms,
|
||||
& ipol, nqpts, nbf_ao, mbf_ao, grad,
|
||||
& Fmat, Pmat, Bmat, Dmat,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), ibf_ao, iniz, ifin,
|
||||
& dbl_mb(k_drhonuc), dbl_mb(k_ddelrhonuc),
|
||||
& dbl_mb(irscr), rchi_atom,
|
||||
& rdelchi_atom, rdens_atom, cetobfr)
|
||||
#endif
|
||||
endif
|
||||
c
|
||||
if (.not.ma_pop_stack(lrscr))
|
||||
& call errquit('xc_quadv0_a: cannot pop stack',0)
|
||||
c
|
||||
if (.not.ma_pop_stack(lbp))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
c Evaluate the integrated charge density.
|
||||
c
|
||||
nqtot = nqtot + nqpts
|
||||
|
||||
rho_int = rho_int +
|
||||
& ddot(nqpts, qwght, 1, rho, 1)
|
||||
|
||||
call xc_eval_fnl(rho, delrho, Amat, Amat2, Cmat, Cmat2,
|
||||
& nqpts, Ex, Ec, qwght, GRAD, ldew, func,
|
||||
& do_hess .or. do_cpks_r)
|
||||
c
|
||||
c
|
||||
c Exchange-correlation contribution
|
||||
c
|
||||
c exchange-correlation contribution, with or without gradient
|
||||
c terms
|
||||
c
|
||||
if (do_force) then
|
||||
c
|
||||
call dftg_quadv0_b(ictr_buf, ldew, AO_bas_han, ipol,
|
||||
& grad, nbf_ao, mbf_ao, natoms, nqpts,
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), Amat, Cmat,
|
||||
& ibf_ao, iniz, ifin,
|
||||
& ppp, ff,
|
||||
& zzz, hh, iga_dens, force,
|
||||
& oactive, nactive)
|
||||
c
|
||||
else if (do_hess) then
|
||||
c
|
||||
if (.not.MA_push_get(mt_dbl, 15*nqpts, 'd2scr',
|
||||
& l_d2scr, k_d2scr)) call errquit
|
||||
& ('dftg_quadv0: not enough stack for d2scr',0)
|
||||
c
|
||||
c Note that k_d2scr is passed in three times, since it is also
|
||||
c used for the scratch arrays ff and s. This is ok since by the
|
||||
c time ff and s are used, the routine is done with the array scr.
|
||||
c -- what about using Bmat? !!!
|
||||
call xc_d2expl(tol_rho, dbl_mb(k_d2scr), Amat, Amat2,
|
||||
& Bmat, Cmat, Cmat2, Dmat, Fmat, Pmat,
|
||||
& dbl_mb(k_d2scr), dbl_mb(k_d2scr+3*nqpts),
|
||||
& dbl_mb(ichi_ao), dbl_mb(idelchi_ao),
|
||||
& dbl_mb(iheschi_ao), dbl_mb(id3chi_ao),
|
||||
& curatoms, ncuratoms, ipol, nqpts,
|
||||
& nbf_ao, mbf_ao, grad, AO_bas_han, natoms, iniz,
|
||||
& ifin, dbl_mb(k_drhonuc), dbl_mb(k_ddelrhonuc),
|
||||
& delrho, iga_dens, hess, ibf_ao,
|
||||
& rchi_atom, rdelchi_atom,
|
||||
& rdens_atom, cetobfr)
|
||||
if (.not.ma_pop_stack(l_d2scr))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
else if (do_cpks_r) then
|
||||
c
|
||||
c The CPKS RHS is conveniently done in two parts: the part
|
||||
c involving the derivs of the density looks just like the CPKS
|
||||
c LHS, so that code is called for that part, while the part
|
||||
c involving only the derivs of the basis functions is calculated
|
||||
c with a separate routine
|
||||
c
|
||||
c First part: LHS-like contribution
|
||||
c
|
||||
c Combine derivs of functional with deriv densities to get coeffs
|
||||
c of bfn pairs in CPKS LHS - these are written back into rho and
|
||||
c delrho
|
||||
c
|
||||
call xc_cpks_coeff(dbl_mb(k_drhonuc),
|
||||
& dbl_mb(k_ddelrhonuc), Amat2, Cmat2, Cmat, delrho,
|
||||
& 3*ncuratoms, ipol, nqpts, GRAD)
|
||||
#if 1
|
||||
c write(*,*)'cpks first part'
|
||||
call xc_tabcd_cpks_lhs(tol_rho, Fmat, Pmat,
|
||||
& dbl_mb(k_drhonuc), Bmat, dbl_mb(k_ddelrhonuc),
|
||||
& Dmat, Dbl_MB(ichi_ao), Dbl_MB(idelchi_ao),
|
||||
& 3*ncuratoms, ipol, nqpts, mbf_ao,
|
||||
& nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2, GRAD,
|
||||
& AO_bas_han, natoms, iniz, ifin, g_rhs,
|
||||
& ibf_ao, rchi_atom,
|
||||
& rdelchi_atom, rdens_atom, cetobfr)
|
||||
c write(*,*)'Cmat after lhs'
|
||||
c call output(cmat,1,nq,1,6,nq,6,1)
|
||||
#endif
|
||||
#if 1
|
||||
c
|
||||
c Second part: Basis function derivative contribution
|
||||
c
|
||||
c write(*,*)'cpks second part'
|
||||
c
|
||||
c Transform Cmat to delrho form for use in xc_cpks_tabcd_rhs
|
||||
c
|
||||
if (grad) then
|
||||
call transform_Cmat(delrho, Cmat, ipol, nqpts)
|
||||
endif
|
||||
call xc_tabcd_cpks_rhs(tol_rho, Fmat, Pmat, Amat,
|
||||
& Bmat, Cmat, Dmat, Dbl_MB(ichi_ao),
|
||||
& Dbl_MB(idelchi_ao), dbl_mb(iheschi_ao),
|
||||
& curatoms,
|
||||
& ncuratoms, ipol, nqpts, mbf_ao,
|
||||
& nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2, GRAD,
|
||||
& AO_bas_han, natoms, iniz, ifin, g_rhs,
|
||||
& ibf_ao, rchi_atom,
|
||||
& rdelchi_atom, rdens_atom, cetobfr)
|
||||
#endif
|
||||
c
|
||||
endif
|
||||
c
|
||||
if (.not.ma_pop_stack(ld3chi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lheschi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lgradchi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
if (.not.ma_pop_stack(lchi_ao))
|
||||
& call errquit('dftg_quadv0a: cannot pop stack',0)
|
||||
c
|
||||
if (do_hess .or. do_cpks_r) then
|
||||
if (grad) then
|
||||
if (.not.ma_pop_stack(l_ddelrhonuc))
|
||||
& call errquit('xc_quadv0_a: cannot pop ddelrhonuc',
|
||||
& 0)
|
||||
endif
|
||||
if (.not.ma_pop_stack(l_drhonuc))
|
||||
& call errquit('xc_quadv0_a: cannot pop drhonuc',0)
|
||||
endif
|
||||
c
|
||||
c
|
||||
cedo no weight derivatives for the moment
|
||||
cedo if (do_force.and.ldew)then
|
||||
cedo write(0,*) ictr_buf,natoms,nqpts
|
||||
cedo do i=1,natoms
|
||||
cedo iandex(ii)=ii
|
||||
cedo enddo
|
||||
cedo call dftg_wderiv(ictr_buf, natoms, nqpts, nqpts,
|
||||
cedo . iandex,
|
||||
cedo & qcenters, func, d_qwght,
|
||||
cedo & force, oactive, nactive)
|
||||
cedo endif
|
||||
c
|
||||
c
|
||||
return
|
||||
c
|
||||
end
|
||||
Loading…
Add table
Add a link
Reference in a new issue