comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-29 06:35:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

HvD: Added the ROKS analytical gradients capability.

Browse source 
This capability essentially exploits the fact that ROKS can be considered a
variant of UKS with additional constraints that the alpha and beta orbitals be
the same. 
Therefore at the top layer the implementation imposes the constraints of ROKS.
Hence we have a single set of orbitals and a single Fock matrix, etc. When it
comes, however, to evaluating the Fock matrix the implementation switches to
using a UKS approach. Hence we have alpha and beta densities (or closed shell
and spin densities depending on how UKS is implemented), and we construct 
alpha and beta Fock matrices. The main difference between the UKS and ROKS
gradient is the Lagrangian. For the ROKS formalism the Lagrangian can be copied
from the ROHF formalism.
The implementation requires nothing more than carefully threading the
appropriate arrangements through the code. The main issue being that in ROKS
ipol=2 whereas in ROHF it is ipol=1. This is related to ROKS being based on UKS
whereas in the ROHF this is handled differently. Maybe this can ultimately be
addressed by handling ROKS and ROHF in a unified way. Obviously in such a
situation the ROKS implementation imposes the most stringent constraints. This
is due to the non-linearity of the density functionals which means that ROKS
cannot be cast in terms of closed shell and open shell densities, but only in
terms of alpha and beta densities. This distinction does not exist at the ROHF
level (because in ROHF the Fock matrices are linear in terms of the density
matrix).
This commit is contained in:
Huub Van Dam 2011-10-08 05:49:59 +00:00
parent 92772c41f8
commit 721ddefac5
16 changed files with 2773 additions and 307 deletions
Download patch file Download diff file Expand all files Collapse all files

2
QA/doqmtests.mpi
View file

@ -31,7 +31,7 @@ endif
./runtests.mpi.unix procs $np dplot dft_meta dft_mpwb1khf dft_m05nh2ch3 prop_uhf_h2o
./runtests.mpi.unix procs $np et_zn_dimer vectors_rotate sad_ch3hf
./runtests.mpi.unix procs $np h2o_diag_to_cg h2o_cg_to_diag
./runtests.mpi.unix procs $np h2o_diag_opt h2o_cg_opt h2o_diag_opt_ub3lyp h2o_cg_opt_ub3lyp
./runtests.mpi.unix procs $np h2o_diag_opt h2o_cg_opt h2o_diag_opt_ub3lyp h2o_cg_opt_ub3lyp h2o_cg_opt_rob3lyp
./runtests.mpi.unix procs $np h2o_diag_to_cg_ub3lyp h2o_cg_to_diag_ub3lyp
#
#--- small tests that should fail!

28
QA/tests/h2o_cg_opt_rob3lyp/h2o_cg_opt_rob3lyp.nw Normal file
View file

@ -0,0 +1,28 @@
echo
start h2o_cg_opt_rob3lyp_dat
memory global 40 mb stack 23 mb heap 5 mb
geometry units au
O 0 0 0
H 0 1.430 -1.107
H 0 -1.430 -1.107
end
charge 1
basis
O library 6-31g*
H library 6-31g*
end
dft
cgmin
mult 2
rodft
xc slater
print "forces" "force components"
convergence nr 0.0
end
task dft optimize

2346
QA/tests/h2o_cg_opt_rob3lyp/h2o_cg_opt_rob3lyp.out Normal file
View file

File diff suppressed because it is too large Load diff

4
src/ddscf/fock_xc.F
View file

@ -210,9 +210,9 @@ c -----------------------
if(ndens.eq.2) then
blo(1) = 2*nfock_xc+1+1
bhi(1) = 2*nfock_xc+1+1
int_mb(igxcd+2*nfock_xc+1) =
int_mb(igxcd+2*nfock_xc+1) =
= ga_create_atom_blocked(geom,ao_bas_han,'gdens+2')
call ga_copy(g_dens_scf(2),int_mb(igxcd+2*nfock_xc+1))
call ga_copy(g_dens_scf(2),int_mb(igxcd+2*nfock_xc+1))
endif
else
g_dens_scf(1) = int_mb(igxcd)

8
src/gradients/grad2.F
View file

@ -277,10 +277,10 @@ c
$ d_jl3, d_il3, d_jk3,
$ d_ij4, d_kl4, d_ik4,
$ d_jl4, d_il4, d_jk4,
$ d_ik5, d_jl5, d_il5, d_jk5,
$ d_ik6, d_jl6, d_il6, d_jk6,
$ d_ik7, d_jl7, d_il7, d_jk7,
$ d_ik8, d_jl8, d_il8, d_jk8,
$ d_ik5, d_jl5, d_il5, d_jk5,
$ d_ik6, d_jl6, d_il6, d_jk6,
$ d_ik7, d_jl7, d_il7, d_jk7,
$ d_ik8, d_jl8, d_il8, d_jk8,
$ blen,
$ iilo, jjlo, kklo, lllo,
$ iihi, jjhi, kkhi, llhi,

502
src/gradients/grad_dens.F
View file

@ -1,6 +1,6 @@
subroutine grad_dens (geom, basis, g_dens, g_wdens, g_eigen_diag,
$ evals, occ, ndens, nbf, nbf2, nopen, nclosed, scftype,
$ movecs, omp2, pdm2, coeff,lsmear,osep)
$ movecs, omp2, odft, pdm2, coeff,lsmear,osep)
c $Id$
@ -26,6 +26,7 @@ C movecs (and Lagrangian in ROHF case) on disk
integer l_buf, k_buf
logical omp2
logical odft ! [input] doing DFT (ROKS needs to know)
logical osep ! [input] get seperate alpha and beta density
! for UHF
@ -67,45 +68,45 @@ C get MO vectors from file (taken from Roberts scf_movecs_read)
c
if (scftype .eq. 'RHF') then
C energy-weighted density
C create eigenvalue diagonal matrix
c
C energy-weighted density
c
C create eigenvalue diagonal matrix
nocc = nclosed
if(lsmear) then
if(lsmear) then
c
c smearing, therefore check new nocs
c smearing, therefore check new nocs
c
if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
c
g_tmp = ga_create_atom_blocked(geom, basis,
& 'frac vecs')
call ga_zero(g_tmp)
do i = ga_nodeid()+1, nbf, ga_nnodes()
if(occ(i).gt.1.d-9) then
call get_col(g_vecs, nbf, i, DBL_MB(itmpm))
call dscal(nbf, occ(i)*0.5d0, DBL_MB(itmpm), 1)
call put_col(g_tmp, nbf, i, DBL_MB(itmpm))
g_tmp = ga_create_atom_blocked(geom, basis,
& 'frac vecs')
call ga_zero(g_tmp)
do i = ga_nodeid()+1, nbf, ga_nnodes()
if(occ(i).gt.1.d-9) then
call get_col(g_vecs, nbf, i, DBL_MB(itmpm))
call dscal(nbf, occ(i)*0.5d0, DBL_MB(itmpm), 1)
call put_col(g_tmp, nbf, i, DBL_MB(itmpm))
endif
enddo
do i=1,nbf
if(occ(i).gt.1.d-9) nocc=i
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
if(nocc.gt.nclosed) then
if(.not.ga_destroy ( g_eigen_diag )) call errquit(
( 'grad_dens: cannot destroy eigendiag',11, GA_ERR)
if (.not. ga_create(mt_dbl, nocc, nocc,'eigen_diag',
, 1, nocc, g_eigen_diag))
, call errquit('gradients: ga diags?',nocc*nocc, GA_ERR)
call ga_zero(g_eigen_diag)
endif
enddo
do i=1,nbf
if(occ(i).gt.1.d-9) nocc=i
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
if(nocc.gt.nclosed) then
if(.not.ga_destroy ( g_eigen_diag )) call errquit(
( 'grad_dens: cannot destroy eigendiag',11, GA_ERR)
if (.not. ga_create(mt_dbl, nocc, nocc,'eigen_diag',
, 1, nocc, g_eigen_diag))
, call errquit('gradients: ga diags?',nocc*nocc, GA_ERR)
call ga_zero(g_eigen_diag)
endif
nclosed=nocc
else
g_tmp=g_vecs
endif
nclosed=nocc
else
g_tmp=g_vecs
endif
#if 0
if (ga_nodeid() .eq. 0) then
status = MA_alloc_get(MT_INT, nocc, 'indices', l_ind, k_ind)
@ -120,10 +121,11 @@ c
status = MA_free_heap ( l_ind2 )
end if
#else
call grad_diagb(g_eigen_diag, evals, nocc)
call grad_diagb(g_eigen_diag, evals, nocc)
#endif
C multiply to weighted density
c
C multiply to weighted density
c
call ga_matmul_patch ('n','n', two, zero,
$ g_tmp, 1, nbf, 1, nocc,
$ g_eigen_diag, 1, nocc, 1, nocc,
@ -134,26 +136,32 @@ C multiply to weighted density
$ g_wdens, 1, nbf, 1, nbf )
if ( omp2 ) then
C read 2nd order correction to weighted density from disk
C use g_dens(3) as buffer
c
C read 2nd order correction to weighted density from disk
C use g_dens(3) as buffer
c
call util_file_name ('w_a_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(3) )
call ga_add ( one, g_wdens, mtwo, g_dens(3), g_wdens )
end if
C density matrix
c
C density matrix
c
call ga_matmul_patch('n', 't', 2.0d0, 0.0d0,
$ g_tmp, 1, nbf, 1, nocc,
$ g_vecs, 1, nocc, 1, nbf,
$ g_dens(1), 1, nbf, 1, nbf)
call ga_symmetrize(g_dens(1))
if ( omp2 ) then
C g_dens(1) = D(HF+(2)) for 1-el. contribution
C g_dens(2) = D(HF) for 2-el. contr.
C g_dens(3) = D(HF+2(2)) "
c
C g_dens(1) = D(HF+(2)) for 1-el. contribution
C g_dens(2) = D(HF) for 2-el. contr.
C g_dens(3) = D(HF+2(2)) "
call ga_copy ( g_dens(1), g_dens(2) )
C read and add MP2 contribution
c
C read and add MP2 contribution
c
call util_file_name ('p_a_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(3) )
call ga_scale ( g_dens(3), two )
@ -164,15 +172,16 @@ C read and add MP2 contribution
else if (scftype .eq. 'ROHF') then
nocc = nopen + nclosed
C read in Lagrangian matrix ('eigenvalue matrix' in ROHF case)
C read in Lagrangian matrix ('eigenvalue matrix' in ROHF case)
call util_file_name('lagr', .true.,.false., fname )
status = file_read_ga( fname, g_eigen_diag )
if (.not.status) then
call errquit('grad_dens: could not read Lagrangian', 110,
& GA_ERR)
end if
C multiply to weighted density
c
C multiply to weighted density
c
call ga_matmul_patch ('n','n', two, zero,
$ g_vecs, 1, nbf, 1, nocc,
$ g_eigen_diag, 1, nocc, 1, nocc,
@ -181,216 +190,257 @@ C multiply to weighted density
$ g_dens(1), 1, nbf, 1, nocc,
$ g_vecs, 1, nocc, 1, nbf,
$ g_wdens, 1, nbf, 1, nbf )
c
if (.not.odft) then
c
C density matrices:
C g_dens(1): D(closed) + 1/2 D(open)
C g_dens(2): D(closed)
C g_dens(3): D(open)
C density contains a factor 2 compared to most literature
C (same as in ddscf (I think...))
c
call ga_matmul_patch('n', 't', 2.0d0, 0.0d0,
$ g_vecs, 1, nbf, 1, nclosed,
$ g_vecs, 1, nclosed, 1, nbf,
$ g_dens(1), 1, nbf, 1, nbf )
C density matrices:
C g_dens(1): D(closed) + 1/2 D(open)
C g_dens(2): D(closed)
C g_dens(3): D(open)
C density contains a factor 2 compared to most literature
C (same as in ddscf (I think...))
call ga_matmul_patch('n', 't', 2.0d0, 0.0d0,
$ g_vecs, 1, nbf, nclosed+1, nocc,
$ g_vecs, nclosed+1, nocc, 1, nbf,
$ g_dens(3), 1, nbf, 1, nbf )
call ga_matmul_patch('n', 't', 2.0d0, 0.0d0,
$ g_vecs, 1, nbf, 1, nclosed,
$ g_vecs, 1, nclosed, 1, nbf,
$ g_dens(1), 1, nbf, 1, nbf )
call ga_copy ( g_dens(1), g_dens(2) )
call ga_dadd (one, g_dens(1), 0.5D0, g_dens(3), g_dens(1))
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
call ga_symmetrize(g_dens(3))
c
else
c
C density matrices:
C g_dens(1): D(alpha)
C g_dens(2): D(beta)
c if (.not.osep) then:
C g_dens(1): D(alpha)+D(beta)
C g_dens(2): D(alpha)-D(beta)
c
call ga_matmul_patch('n', 't', 1.0d0, 0.0d0,
$ g_vecs, 1, nbf, 1, nocc,
$ g_vecs, 1, nocc, 1, nbf,
$ g_dens(1), 1, nbf, 1, nbf )
call ga_matmul_patch('n', 't', 2.0d0, 0.0d0,
$ g_vecs, 1, nbf, nclosed+1, nocc,
$ g_vecs, nclosed+1, nocc, 1, nbf,
$ g_dens(3), 1, nbf, 1, nbf )
call ga_matmul_patch('n', 't', 1.0d0, 0.0d0,
$ g_vecs, 1, nbf, 1, nclosed,
$ g_vecs, 1, nclosed, 1, nbf,
$ g_dens(2), 1, nbf, 1, nbf )
call ga_copy ( g_dens(1), g_dens(2) )
call ga_dadd (one, g_dens(1), 0.5D0, g_dens(3), g_dens(1))
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
call ga_symmetrize(g_dens(3))
if (.not. osep) then
call ga_dadd ( one, g_dens(1), one, g_dens(2), g_dens(1) )
call ga_dadd ( one, g_dens(1), mtwo, g_dens(2), g_dens(2) )
endif
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
c
endif
c
else if (scftype .eq. 'UHF') then
C energy-weighted density
C create eigenvalue diagonal matrix
c
C energy-weighted density
c
C create eigenvalue diagonal matrix
nocca = nclosed + nopen ! == nalpha
noccb = nclosed ! nbeta
if(lsmear) then
if(lsmear) then
c
c smearing, therefore check new nocs
c smearing, therefore check new nocs
c
if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
c
g_tmp = ga_create_atom_blocked(geom, basis,
& 'frac vecs')
call ga_zero(g_tmp)
do i = ga_nodeid()+1, nbf, ga_nnodes()
if(occ(i).gt.1.d-9) then
call get_col(g_vecs, nbf, i, DBL_MB(itmpm))
call dscal(nbf, occ(i), DBL_MB(itmpm), 1)
call put_col(g_tmp, nbf, i, DBL_MB(itmpm))
endif
enddo
do i=1,nbf
if(occ(i).gt.1.d-9) nocca=i
enddo
do i=nbf+1,nbf*2
if(occ(i).gt.1.d-9) noccb=i-nbf
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
g_tmp = ga_create_atom_blocked(geom, basis,
& 'frac vecs')
call ga_zero(g_tmp)
do i = ga_nodeid()+1, nbf, ga_nnodes()
if(occ(i).gt.1.d-9) then
call get_col(g_vecs, nbf, i, DBL_MB(itmpm))
call dscal(nbf, occ(i), DBL_MB(itmpm), 1)
call put_col(g_tmp, nbf, i, DBL_MB(itmpm))
endif
enddo
do i=1,nbf
if(occ(i).gt.1.d-9) nocca=i
enddo
do i=nbf+1,nbf*2
if(occ(i).gt.1.d-9) noccb=i-nbf
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
else
g_tmp=g_vecs
endif
if(nocca.gt.nopen+nclosed) then
else
g_tmp=g_vecs
endif
if(nocca.gt.nopen+nclosed) then
if(.not.ga_destroy ( g_eigen_diag )) call errquit(
( 'grad_dens: cannot destroy eigendiag',11, GA_ERR)
& 'grad_dens: cannot destroy eigendiag',11, GA_ERR)
if (.not. ga_create(mt_dbl, nocca, nocca,'eigen_diag',
, 1, nocca, g_eigen_diag))
, call errquit('gradients: ga diags?',nocca*nocca, GA_ERR)
, 1, nocca, g_eigen_diag))
, call errquit('gradients: ga diags?',nocca*nocca, GA_ERR)
call ga_zero(g_eigen_diag)
endif
#if 0
if (ga_nodeid() .eq. 0) then
mxnoc=max(nocca,noccb)
status=MA_alloc_get(MT_INT,mxnoc, 'indices', l_ind, k_ind)
status=MA_alloc_get(MT_INT,mxnoc, 'indices',l_ind2,k_ind2)
do i=1, nocca
int_mb( k_ind+i-1 ) = i
int_mb( k_ind2+i-1 ) = i
end do
call ga_scatter ( g_eigen_diag, evals, int_mb(k_ind),
$ int_mb(k_ind2), nocca)
end if
if (ga_nodeid() .eq. 0) then
mxnoc=max(nocca,noccb)
status=MA_alloc_get(MT_INT,mxnoc, 'indices', l_ind, k_ind)
status=MA_alloc_get(MT_INT,mxnoc, 'indices',l_ind2,k_ind2)
do i=1, nocca
int_mb( k_ind+i-1 ) = i
int_mb( k_ind2+i-1 ) = i
end do
call ga_scatter ( g_eigen_diag, evals, int_mb(k_ind),
$ int_mb(k_ind2), nocca)
end if
#else
call grad_diagb(g_eigen_diag, evals, nocca)
call grad_diagb(g_eigen_diag, evals, nocca)
#endif
C multiply to weighted density
call ga_matmul_patch ('n','n', one, zero,
$ g_tmp, 1, nbf, 1, nocca,
$ g_eigen_diag, 1, nocca, 1, nocca,
$ g_dens(1), 1, nbf, 1, nocca )
call ga_matmul_patch ('n','t', one, zero,
$ g_dens(1), 1, nbf, 1, nocca,
$ g_vecs, 1, nocca, 1, nbf,
$ g_wdens, 1, nbf, 1, nbf )
if(lsmear) then
c
c smearing, therefore check new nocs
C multiply to weighted density
call ga_matmul_patch ('n','n', one, zero,
$ g_tmp, 1, nbf, 1, nocca,
$ g_eigen_diag, 1, nocca, 1, nocca,
$ g_dens(1), 1, nbf, 1, nocca )
call ga_matmul_patch ('n','t', one, zero,
$ g_dens(1), 1, nbf, 1, nocca,
$ g_vecs, 1, nocca, 1, nbf,
$ g_wdens, 1, nbf, 1, nbf )
if(lsmear) then
c
if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
c smearing, therefore check new nocs
c
if (.not. MA_Push_Get(MT_Dbl, nbf, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
c
g_tmp2 = ga_create_atom_blocked(geom, basis,
& 'frac vecs')
& 'frac vecs')
call ga_zero(g_tmp2)
do i = ga_nodeid()+1, nbf, ga_nnodes()
if(occ(i).gt.1.d-9) then
call get_col(g_vecs2, nbf, i, DBL_MB(itmpm))
call dscal(nbf, occ(i+nbf), DBL_MB(itmpm), 1)
call put_col(g_tmp2, nbf, i, DBL_MB(itmpm))
endif
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
if(noccb.gt.nclosed) then
if(.not.ga_destroy ( g_eigen_diag )) call errquit(
( 'grad_dens: cannot destroy eigendiag',11, GA_ERR)
if (.not. ga_create(mt_dbl, noccb, noccb,'eigen_diag',
, 1, noccb, g_eigen_diag))
, call errquit('gradients: ga diags?',noccb*noccb, GA_ERR)
call ga_zero(g_eigen_diag)
endif
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
if(noccb.gt.nclosed) then
if(.not.ga_destroy ( g_eigen_diag )) call errquit(
& 'grad_dens: cannot destroy eigendiag',11, GA_ERR)
if (.not. ga_create(mt_dbl, noccb, noccb,'eigen_diag',
& 1, noccb, g_eigen_diag))
& call errquit('gradients: ga diags?',noccb*noccb,
& GA_ERR)
call ga_zero(g_eigen_diag)
endif
else
g_tmp2=g_vecs2
endif
else
g_tmp2=g_vecs2
endif
#if 0
if (ga_nodeid() .eq. 0) then
do i=1, noccb
int_mb( k_ind+i-1 ) = i
int_mb( k_ind2+i-1 ) = i
end do
call ga_scatter ( g_eigen_diag, evals(nbf+1), int_mb(k_ind),
$ int_mb(k_ind2), noccb )
status = MA_free_heap ( l_ind )
status = MA_free_heap ( l_ind2 )
end if
if (ga_nodeid() .eq. 0) then
do i=1, noccb
int_mb( k_ind+i-1 ) = i
int_mb( k_ind2+i-1 ) = i
end do
call ga_scatter ( g_eigen_diag, evals(nbf+1), int_mb(k_ind),
$ int_mb(k_ind2), noccb )
status = MA_free_heap ( l_ind )
status = MA_free_heap ( l_ind2 )
end if
#else
call grad_diagb(g_eigen_diag, evals(nbf+1), noccb)
call grad_diagb(g_eigen_diag, evals(nbf+1), noccb)
#endif
C multiply to weighted density and add to spin up part
call ga_matmul_patch ('n','n', one, zero,
$ g_tmp2, 1, nbf, 1, noccb,
$ g_eigen_diag, 1, noccb, 1, noccb,
$ g_dens(1), 1, nbf, 1, noccb )
call ga_matmul_patch ('n','t', one, one,
$ g_dens(1), 1, nbf, 1, noccb,
$ g_vecs2, 1, noccb, 1, nbf,
$ g_wdens, 1, nbf, 1, nbf )
call ga_symmetrize(g_wdens)
C multiply to weighted density and add to spin up part
call ga_matmul_patch ('n','n', one, zero,
$ g_tmp2, 1, nbf, 1, noccb,
$ g_eigen_diag, 1, noccb, 1, noccb,
$ g_dens(1), 1, nbf, 1, noccb )
call ga_matmul_patch ('n','t', one, one,
$ g_dens(1), 1, nbf, 1, noccb,
$ g_vecs2, 1, noccb, 1, nbf,
$ g_wdens, 1, nbf, 1, nbf )
call ga_symmetrize(g_wdens)
if ( omp2 ) then
C read 2nd order correction to weighted density from disk
C use g_dens(3) as buffer
call util_file_name ('w_a_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(3) )
call ga_add ( one, g_wdens, mone, g_dens(3), g_wdens )
call util_file_name ('w_b_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(3) )
call ga_add ( one, g_wdens, mone, g_dens(3), g_wdens )
end if
if ( omp2 ) then
c
C read 2nd order correction to weighted density from disk
C use g_dens(3) as buffer
c
call util_file_name ('w_a_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(3) )
call ga_add ( one, g_wdens, mone, g_dens(3), g_wdens )
call util_file_name ('w_b_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(3) )
call ga_add ( one, g_wdens, mone, g_dens(3), g_wdens )
end if
c
C density matrix
C D_up and D_down (or and b)
call ga_matmul_patch('n', 't', 1.0d0, 0.0d0,
$ g_tmp, 1, nbf, 1, nocca,
$ g_vecs, 1, nocca, 1, nbf,
$ g_dens(1), 1, nbf, 1, nbf)
call ga_matmul_patch('n', 't', 1.0d0, 0.0d0,
$ g_tmp2, 1, nbf, 1, noccb,
$ g_vecs2, 1, noccb, 1, nbf,
$ g_vecs, 1, nbf, 1, nbf)
C density matrix
C D_up and D_down (or and b)
call ga_matmul_patch('n', 't', 1.0d0, 0.0d0,
$ g_tmp, 1, nbf, 1, nocca,
$ g_vecs, 1, nocca, 1, nbf,
$ g_dens(1), 1, nbf, 1, nbf)
call ga_matmul_patch('n', 't', 1.0d0, 0.0d0,
$ g_tmp2, 1, nbf, 1, noccb,
$ g_vecs2, 1, noccb, 1, nbf,
$ g_vecs, 1, nbf, 1, nbf)
if ( .not. omp2 ) then
C D+ and D- (Formalism from Dupuis & King, JCP 68(9), p4000
if (.not. osep) then
if ( .not. omp2 ) then
c
C D+ and D- (Formalism from Dupuis & King, JCP 68(9), p4000
c
if (.not. osep) then
call ga_dadd ( one, g_dens(1), mone, g_vecs, g_dens(2) )
call ga_dadd ( one, g_dens(1), one, g_vecs, g_dens(1) )
else
call ga_copy (g_vecs,g_dens(2))
endif
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
else ! UMP2
C g_dens(1) = D(a+b) (HF+(2)) for 1-el. contribution
C g_dens(2) = D(a) (HF) for 2-el. contr.
C g_dens(3) = D(b) (HF) "
C g_dens(4) = D(a) (HF+2(2)) "
C g_dens(5) = D(b) (HF+2(2)) "
else
call ga_copy (g_vecs,g_dens(2))
endif
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
else ! UMP2
c
C g_dens(1) = D(a+b) (HF+(2)) for 1-el. contribution
C g_dens(2) = D(a) (HF) for 2-el. contr.
C g_dens(3) = D(b) (HF) "
C g_dens(4) = D(a) (HF+2(2)) "
C g_dens(5) = D(b) (HF+2(2)) "
c
C move HF matrices to g_dens(2:3)
call ga_copy ( g_dens(1), g_dens(2) )
call ga_copy ( g_vecs, g_dens(3) )
C move HF matrices to g_dens(2:3)
call ga_copy ( g_dens(1), g_dens(2) )
call ga_copy ( g_vecs, g_dens(3) )
call ga_dadd ( one, g_dens(1), one, g_dens(3), g_dens(1) )
C read and add MP2 contribution for spin up
call util_file_name ('p_a_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(4) )
call ga_dadd ( one, g_dens(4), one, g_dens(1), g_dens(1) )
call ga_dadd ( two, g_dens(4), one, g_dens(2), g_dens(4) )
C spin down
call util_file_name ('p_b_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(5) )
call ga_dadd ( one, g_dens(5), one, g_dens(1), g_dens(1) )
call ga_dadd ( two, g_dens(5), one, g_dens(3), g_dens(5) )
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
call ga_symmetrize(g_dens(3))
call ga_symmetrize(g_dens(4))
call ga_symmetrize(g_dens(5))
end if ! UMP2
call ga_dadd ( one, g_dens(1), one, g_dens(3), g_dens(1) )
c
C read and add MP2 contribution for spin up
c
call util_file_name ('p_a_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(4) )
call ga_dadd ( one, g_dens(4), one, g_dens(1), g_dens(1) )
call ga_dadd ( two, g_dens(4), one, g_dens(2), g_dens(4) )
c
C spin down
c
call util_file_name ('p_b_ao', .true., .true., fname )
status = file_read_ga ( fname, g_dens(5) )
call ga_dadd ( one, g_dens(5), one, g_dens(1), g_dens(1) )
call ga_dadd ( two, g_dens(5), one, g_dens(3), g_dens(5) )
call ga_symmetrize(g_dens(1))
call ga_symmetrize(g_dens(2))
call ga_symmetrize(g_dens(3))
call ga_symmetrize(g_dens(4))
call ga_symmetrize(g_dens(5))
end if ! UMP2
c
else if (scftype .eq. 'MCSCF') then
c
@ -398,35 +448,35 @@ c
call util_file_name('mcpdm',.true.,.false.,pdmfilename)
if (eaf_open(pdmfilename, eaf_r, pdmfile) .ne. 0)
$ call errquit('mcscf grad: failed opening 2pdm file',0,
& GA_ERR)
& GA_ERR)
offset = 0.0d0
if (.not. ma_push_get(mt_dbl,nbf*nbf,'buffer',l_buf,k_buf))
$ call errquit('mcscf grad: memory problem',nbf**2,
& MA_ERR)
& MA_ERR)
if (eaf_read(pdmfile, offset, dbl_mb(k_buf), 8*nbf**2).ne.0)
$ call errquit('mcscf grad: failed reading lagrangian',0,
& DISK_ERR)
& DISK_ERR)
call ga_put(g_wdens, 1, nbf, 1, nbf, dbl_mb(k_buf), nbf)
offset = offset + 8*nbf**2
if (eaf_read(pdmfile, offset, dbl_mb(k_buf), 8*nbf**2).ne.0)
$ call errquit('mcscf grad: failed reading 1-density',0,
& DISK_ERR)
& DISK_ERR)
call ga_put(g_dens(1), 1, nbf, 1, nbf, dbl_mb(k_buf), nbf)
offset = offset + 8*nbf**2
if (eaf_read(pdmfile, offset, pdm2, 8*nopen**4).ne.0)
$ call errquit('mcscf grad: failed reading 2pdm ',0,
& DISK_ERR)
& DISK_ERR)
call ga_get(g_vecs, 1, nbf, nclosed+1, nclosed+nopen,
$ coeff, nbf)
if (eaf_close(pdmfile).ne.0)
$ call errquit('mcscf: failed closing 2pdm file',0,
& DISK_ERR)
& DISK_ERR)
* call MA_summarize_allocated_blocks()

37
src/gradients/grad_force.F
View file

@ -132,7 +132,7 @@ c
character*16 tag
character*32 theory
character*8 scftype
character*8 scftype, lscftype ! lscftype is just for calling grad2
character*32 rtdb_string
c
@ -215,7 +215,8 @@ C get SCF MO vectors for density
if (.not. rtdb_cget(rtdb, 'mcscf:input vectors', 1, movecs_in))
$ call errquit('gradients: MCSCF MO vectors not defined',0,
& RTDB_ERR)
scftype = 'mcscf'
scftype = 'mcscf'
lscftype = 'mcscf'
else if (theory .eq. 'dft') then
odft = .true.
if (.not. rtdb_cget(rtdb, 'dft:input vectors', 1, movecs_in))
@ -228,11 +229,11 @@ C get SCF MO vectors for density
& RTDB_ERR)
if (.not. rtdb_get(rtdb, 'dft:ipol', mt_int, 1, ipol))
$ ipol = 1
if (ipol .eq. 1) then
scftype = 'rhf'
else
scftype = 'uhf'
endif
if (.not. rtdb_cget(rtdb, 'dft:scftype', 1, scftype))
& call errquit('gradients: scftype not defined',0,
& RTDB_ERR)
lscftype = scftype
if (scftype.eq.'ROHF') lscftype = 'UHF'
jfac = 1.0d0
if (ocdfit) jfac = 0.0d0
kfac = xfac(1)
@ -240,6 +241,7 @@ C get SCF MO vectors for density
if (.not. rtdb_cget(rtdb, 'scf:scftype', 1, scftype))
$ call errquit('gradients: scftype not defined',0,
& RTDB_ERR)
lscftype = scftype
if (.not. rtdb_cget(rtdb, 'scf:input vectors', 1, movecs_in))
$ call errquit('gradients: SCF MO vectors not defined',0,
& RTDB_ERR)
@ -267,12 +269,17 @@ C scftype: MCSCF, RHF, ROHF or UHF
if (.not. rtdb_get(rtdb, 'dft:noc', mt_int, 2, noc))
& call errquit('gradients: rtdb_get of noc failed', 0,
& RTDB_ERR)
if (ipol .eq. 1) then
if (scftype.eq.'RHF') then
nclosed = noc(1)
nopen = 0
else
elseif (scftype.eq.'ROHF') then
nclosed = min(noc(1),noc(2))
nopen = abs(noc(1) - noc(2))
elseif (scftype.eq.'UHF') then
nclosed = min(noc(1),noc(2))
nopen = abs(noc(1) - noc(2))
else
call errquit("grad_force: invalid scftype",0,UERR)
endif
else
if (.not. rtdb_get(rtdb, 'scf:nopen', mt_int, 1, nopen))
@ -326,6 +333,10 @@ C # of eigenvalues and density matrices
else ! ROHF
nbf2 = nbf
ndens = 3
if (odft) then
nbf2 = 2*nbf
ndens = 2
endif
end if
end if
@ -463,7 +474,7 @@ c
& frac_occ )) frac_occ=.false.
call grad_dens(geom, basis, g_dens, g_wdens, g_eigen_diag,
$ dbl_mb(k_evals), dbl_mb(k_occ), ndens, nbf, nbf2,
$ nopen, nclosed, scftype, movecs_in, omp2,
$ nopen, nclosed, scftype, movecs_in, omp2, odft,
$ dbl_mb(k_pdm2), dbl_mb(k_coeff),frac_occ,.false.)
C free temporary arrays
@ -853,7 +864,7 @@ c
$ lscratch, dbl_mb(k_scr), lbuf/12, dbl_mb(k_buf),
$ int_mb(k_labels), maxq, int_mb(k_list), dbl_mb(k_q4),
$ tol2e, nsh,
$ log_mb(k_act), oskel, scftype, omp2, nopen, nbf,
$ log_mb(k_act), oskel, lscftype, omp2, nopen, nbf,
$ dbl_mb(k_pdm2),dbl_mb(k_pdm2a),dbl_mb(k_pdm2b), ! MCSCF
$ dbl_mb(k_pdm2c),dbl_mb(k_pdm2d), dbl_mb(k_coeff), ! MCSCF
$ nshblocks,
@ -887,7 +898,7 @@ c
$ lscratch, dbl_mb(k_scr), lbuf/12, dbl_mb(k_buf),
$ int_mb(k_labels), maxq, int_mb(k_list), dbl_mb(k_q4),
$ tol2e, nsh,
$ log_mb(k_act), oskel, scftype, omp2, nopen, nbf,
$ log_mb(k_act), oskel, lscftype, omp2, nopen, nbf,
$ dbl_mb(k_pdm2),dbl_mb(k_pdm2a),dbl_mb(k_pdm2b), ! MCSCF
$ dbl_mb(k_pdm2c),dbl_mb(k_pdm2d), dbl_mb(k_coeff), ! MCSCF
$ nshblocks,
@ -920,7 +931,7 @@ c
$ lscratch, dbl_mb(k_scr), lbuf/12, dbl_mb(k_buf),
$ int_mb(k_labels), maxq, int_mb(k_list), dbl_mb(k_q4),
$ tol2e, nsh,
$ log_mb(k_act), oskel, scftype, omp2, nopen, nbf,
$ log_mb(k_act), oskel, lscftype, omp2, nopen, nbf,
$ dbl_mb(k_pdm2),dbl_mb(k_pdm2a),dbl_mb(k_pdm2b), ! MCSCF
$ dbl_mb(k_pdm2c),dbl_mb(k_pdm2d), dbl_mb(k_coeff), ! MCSCF
$ nshblocks,

6
src/hessian/analytic/hess_init.F
View file

@ -27,7 +27,7 @@ c
c integer rtdb, mult, g_vecs(2)
integer rtdb, mult
integer dftnoc(2), nbsize
logical oprint, omp2
logical oprint, omp2, odft
logical cdfit, frac_occ
integer nbf_cd, nbf2, mxnocc
integer g_eigen_diag
@ -53,6 +53,7 @@ c
c
c get geometry and basis set handles and associated useful information
c
odft = .false.
if (oprint) write(LuOut,*) 'HESS: set up geometry and basis'
if (.not. geom_create(geom, 'geometry'))
* call errquit('hess_init: geom_create problem', 555, GEOM_ERR)
@ -181,6 +182,7 @@ c nopen = nopen + 1
kfac(2) = 1.0d0
endif
elseif (theory .eq. 'dft') then
odft = .true.
if (.not. rtdb_get(rtdb, 'dft:mult', MT_INT, 1,mult))
* mult = 1
nopen = mult - 1
@ -324,7 +326,7 @@ c
if (oprint) write(LuOut,*) 'HESS: getting densities'
call grad_dens(geom, basis, g_dens, g_wdens, g_eigen_diag,
$ dbl_mb(k_evals), dbl_mb(k_occ), ndens, nbf, nbf2,
$ nopen, nclosed, scftype, movecs_in, omp2,
$ nopen, nclosed, scftype, movecs_in, omp2, odft,
$ dbl_mb(k_pdm2), dbl_mb(k_coeff),frac_occ,.true.)
c
c free temporary arrays

86
src/nwdft/dftgrad/dft_gradients.F
View file

@ -65,6 +65,7 @@ c integer noc(2)
external grad_norm, grad_max
character*255 title_vecs, basis_vecs
character*20 scftype_vecs
character*80 scftype
logical status,frac_occ
integer me, nproc, max_sh_bf, max_at_bf, nat, max_sh_bfcd,
$ lforce, nactive, i, nbf_vecs, nsets, ispin,
@ -226,18 +227,40 @@ c
$ call errquit('dft_gradient: DFT MO vectors not defined',0,
& RTDB_ERR)
status = movecs_read_header(movecs_in, title_vecs, basis_vecs,
$ scftype_vecs, nbf_vecs, nsets, nmo, 2)
c
if (ipol .ne. nsets)then ! 1 for RHF, 2 for UHF
write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets
call errquit('dft_gradient: ERROR ipol, nsets disagree',2,
& INPUT_ERR)
endif
$ scftype_vecs, nbf_vecs, nsets, nmo, 2)
c
c ipol - number of spin channels (RKS=1, ROKS=2, UKS=2)
c nsets - number of sets of vectors (RKS=1, ROKS=1, UKS=2)
c
if (.not. rtdb_cget(rtdb, 'dft:scftype', 1, scftype))
$ call errquit('dft_gradient: DFT scftype not defined',0,
& RTDB_ERR)
if (scftype.eq.'RHF') then
if (ipol .ne. nsets .or. ipol.ne.1)then
write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets
call errquit('dft_gradient: ERROR ipol, nsets disagree',2,
& INPUT_ERR)
endif
elseif (scftype.eq.'ROHF') then
if (nsets.ne.1.or.ipol.ne.2) then
write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets
call errquit('dft_gradient: ERROR ipol, nsets disagree',2,
& INPUT_ERR)
endif
elseif (scftype.eq.'UHF') then
if (nsets.ne.2.or.ipol.ne.2) then
write (LuOut,*) 'dft_gradient: ERROR ipol, nsets:',ipol,nsets
call errquit('dft_gradient: ERROR ipol, nsets disagree',2,
& INPUT_ERR)
endif
else
call errquit('dft_gradient: illegal scftype',0,UERR)
endif
c
c Should check much more info than just nbf for consistency
c
c
c get mo eigevectors
c get mo eigenvectors
c
if (nbf_ao .ne. nbf_vecs)then
write(LuOut,*)'dft_gradient movecs output = ',movecs_in
@ -247,9 +270,9 @@ c
status = .true.
do ispin = 1, ipol
status = status .and.
$ movecs_read(movecs_in, ispin,
& dbl_mb(k_occ+(ispin-1)*nbf_ao),
$ dbl_mb(k_evals(ispin)), g_vecs(ispin))
$ movecs_read(movecs_in, min(ispin,nsets),
& dbl_mb(k_occ+(ispin-1)*nbf_ao),
$ dbl_mb(k_evals(ispin)), g_vecs(ispin))
enddo
endif
c
@ -261,11 +284,11 @@ c
c
if(frac_occ) then
c
c smearing, therefore check new nocs
c smearing, therefore check new nocs
c
if (.not. MA_Push_Get(MT_Dbl, nbf_ao, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
rhffact = one
if (.not. MA_Push_Get(MT_Dbl, nbf_ao, 'tmpm', ltmpm, itmpm))
& call errquit('dftgforce: failed to alloc tmpm',0, MA_ERR)
rhffact = one
c
do ispin=1,ipol
g_tmp(ispin) = ga_create_atom_blocked(geom, ao_bas_han,
@ -283,8 +306,8 @@ c
if(dbl_mb(ipoint+i).gt.toll) noc(ispin)=i
enddo
enddo
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
if (.not.ma_pop_stack(ltmpm))
& call errquit('dftg_force: cannot pop stack',0, MA_ERR)
else
do ispin=1,ipol
@ -293,12 +316,12 @@ c
endif
do ispin = 1, ipol
c
c dens = vecs*vecs
c dens = vecs*vecs
c
if (odftps) call pstat_on(ps_dgemm)
if (odftps) call pstat_on(ps_dgemm)
call ga_dgemm('n', 't', nbf_ao, nbf_ao, noc(ispin), rhffact,
$ g_tmp(ispin), g_vecs(ispin), 0.0d0, iga_dens(ispin))
if (odftps) call pstat_off(ps_dgemm)
if (odftps) call pstat_off(ps_dgemm)
call ga_symmetrize(iga_dens(ispin))
c
c free temporary arrays
@ -420,10 +443,10 @@ c
$ call errquit('dftg_force: could not alloc svec',0, MA_ERR)
ippp=k_wdens
c
if (odftps) call pstat_on(ps_vcoul)
call xc_rep_init(rtdb, geom, ao_bas_han,iga_dens,iga_dens,
& nbf_ao,ipol,.true.,.true.)
call dftg_cdfit(ao_bas_han, cd_bas_han,
if (odftps) call pstat_on(ps_vcoul)
call xc_rep_init(rtdb, geom, ao_bas_han,iga_dens,iga_dens,
& nbf_ao,ipol,.true.,.true.)
call dftg_cdfit(ao_bas_han, cd_bas_han,
$ nbf_cd, nat, tol2e, dbl_mb(k_scr),
$ lscratch, dbl_mb(k_buf), lbuf,
$ dbl_mb(isvec), dbl_mb(ippp), max_sh_bf,
@ -433,17 +456,18 @@ c
. errquit(' dftggrad: xcrepclose failed ',0, 0)
c
call ga_dgop(msg_grad_2el, dbl_mb(k_frc_2el), lforce, '+')
if (odftps) call pstat_off(ps_vcoul)
if (odftps) call pstat_off(ps_vcoul)
c
if (.not.ma_chop_stack(l_buf))
if (.not.ma_chop_stack(l_buf))
$ call errquit('dft_gradient: cannot chop stack',0, MA_ERR)
c
c restore alpha DM in g_dens(1)
c restore alpha DM in g_dens(1)
c
if (ipol .eq. 2)then
onem = -1.d0
call ga_dadd(one, iga_dens(1), onem, iga_dens(2), iga_dens(1))
endif
if (ipol .eq. 2)then
onem = -1.d0
call ga_dadd(one, iga_dens(1), onem, iga_dens(2),
& iga_dens(1))
endif
endif
c
c get exchange-correlation contribution to the gradient

24
src/nwdft/dftgrad/dftg_getxc.F
View file

@ -61,20 +61,20 @@ c
if (xc_gotxc())then
c
c Compute via direct numerical quadrature.
if(.not.grid_written) then
if(.not.grid_written) then
c
c check if we can get the grid from a file
c check if we can get the grid from a file
c
grid_written=grid_reopen(geom)
endif
if(grid_written) then
call grid_setquad(rtdb)
call grid_packing_info
endif
call dftg_gridv0(rtdb, geom,ao_bas_han,ipol,nbf_ao,
. iga_dens, force, natom, oactive,
& nactive, rdens_atom, cetobfr)
oreopen=.false.
grid_written=grid_reopen(geom)
endif
if(grid_written) then
call grid_setquad(rtdb)
call grid_packing_info
endif
call dftg_gridv0(rtdb, geom,ao_bas_han,ipol,nbf_ao,
. iga_dens, force, natom, oactive,
& nactive, rdens_atom, cetobfr)
oreopen=.false.
c!!! BGJ test !!!
c

2
src/nwdft/dftgrad/dftg_quadv0b.F
View file

@ -1,4 +1,4 @@
Subroutine dftg_quadv0_b(ictr, ldew, AO_bas_han, ipol, GRAD, kske,
subroutine dftg_quadv0_b(ictr, ldew, AO_bas_han, ipol, GRAD, kske,
& nbf_ao, mbf_ao, nat, nq, chi, delchi,
& heschi, Amat, Cmat, Mmat, ccrap,
I ibf, iniz, ifin, cetobfr,

11
src/nwdft/input_dft/dft_inpana.F
View file

@ -137,7 +137,9 @@ c
c odd # of electrons or not a singlet state --> LSD
c
if ((.not.even).or.(mult.ne.1).or.(theory.eq.'sodft').or.
& (scftype.eq.'UHF')) ipol=2
& (scftype.eq.'UHF').or.(scftype.eq.'ROHF')) then
ipol=2
endif
noc(2)=0
c
c Calculate number of occupied orbitals.
@ -336,9 +338,11 @@ c
call util_print_centered
& (LuOut,'General Information',20,.true.)
write(LuOut,9020)
if (ipol.eq.1)then
if (scftype.eq."RHF")then
write(LuOut,9050)
elseif (ipol.eq.2)then
elseif (scftype.eq."ROHF")then
write(LuOut,9052)
elseif (scftype.eq."UHF")then
write(LuOut,9055)
endif
write(LuOut,8150)ncenters
@ -670,6 +674,7 @@ c
9040 format(10x,'Convergence on density requested: ',1Pd8.2)
9045 format(10x,'Convergence on gradient requested: ',1Pd8.2)
9050 format(10x,'Wavefunction type: closed shell.')
9052 format(10x,'Wavefunction type: restricted open shell.')
9055 format(10x,'Wavefunction type: spin polarized.')
9056 format(10x,'Use of symmetry is: ',a3,
& '; symmetry adaption is: ',a3)

2
src/nwdft/scf_dft_cg/dft_cg_info.F
View file

@ -187,7 +187,7 @@ c $ nopen = 0
if (nopen .gt. nelec) call errquit
$ ('scf: nopen > nelec ', nopen, INPUT_ERR)
if (mod(nelec-nopen,2) .ne. 0) call errquit
$ ('scf: no. of closed-shell electrons is not even!',nopen,
$ ('dft: no. of closed-shell electrons is not even!',nopen,
& INPUT_ERR)
c nclosed = (nelec-nopen) / 2
c

7
src/nwdft/scf_dft_cg/dft_roks_fock.F
View file

@ -151,7 +151,12 @@ c cfock = closed shell AO fock
c odens = open shell density
c ocoul = open shell coulomb
c oexch = open shell exchange
c
c
g_a_xc = 0
g_b_xc = 0
g_adens = 0
g_bdens = 0
c
g_cdens = ga_create_atom_blocked(geom, basis, 'rohf: cdens')
g_cfock = ga_create_atom_blocked(geom, basis, 'rohf: cfock')
g_hcore = ga_create_atom_blocked(geom, basis, 'rohf: hcore')

5
src/task/task_gradient.F
View file

@ -289,11 +289,6 @@ c == restricted open-shell DFT ==
if(theory.eq.'dft') then
if (.not. rtdb_get(rtdb,'dft:rodft',mt_log,1,rodft))
1 rodft = .false.
if (rodft) then
numerical = .true. ! do numerical gradients for restricted-open shell
if (ga_nodeid().eq.0) write(LuOut,*)
1 'RODFT - numerical derivatives only'
end if
end if
c
c Currently cosmo does not work with analytic gradients

10
src/util/bgj.F
View file

@ -258,10 +258,10 @@ cDEBUG
status = movecs_read_header(movecs_in, title_vecs, basis_vecs,
& scftype_vecs, nbf_vecs, nsets, nmo, 2)
c
c ipol 1 closed-shell (RKS)
c ipol 2 open-shell (both ROKS and UKS)
c nsets 1 RKS and ROKS
c nsets 2 UKS
c ipol 1 closed-shell (RKS)
c ipol 2 open-shell (both ROKS and UKS)
c nsets 1 RKS and ROKS
c nsets 2 UKS
c
if (ipol .ne. nsets)then
if (ipol.eq.1.and.nsets.eq.2) then
@ -302,7 +302,7 @@ c
& 'bgj_get_scf_dens:beta density')
endif
c
c get mo eigevectors
c get mo eigenvectors
c
if (nbf_ao .ne. nbf_vecs)then
write(6,*)'bgj_get_scf_dens movecs output = ',movecs_in