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

Merge pull request #1263 from edoapra/jmol-fix

Browse source 
Foster-Boys updates
This commit is contained in:
NWChem: Open Source High-Performance Computational Chemistry 2026-03-25 17:34:16 -07:00 • committed by GitHub
commit 179015c5af
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
4 changed files with 232 additions and 190 deletions
Download patch file Download diff file Expand all files Collapse all files

37
src/ddscf/localize.F
View file

@ -227,6 +227,7 @@ c
#include "geom.fh"
#include "bas.fh"
#include "util.fh"
#include "stdio.fh"
c
c Localize the nloc orbitals in iloc() by mixing with each other
c
@ -251,10 +252,10 @@ c locpow (localization sum exponent) should be 2 or 4
gamma_tol = 1d-10
c
if (ga_nodeid() .eq. 0) then
write(6,9001)
write(luout,9001)
9001 format(/10x,' iter Max. delocal Mean delocal Converge'/
& 10x,' ---- ------------ ------------ --------')
call util_flush(6)
call util_flush(luout)
end if
c
dprev = 0.0d0
@ -281,7 +282,7 @@ c
end do
qs = qs + qas**2
end do
* write(6,*) ' ds ', s, 1.0d0/qs
* write(luout,*) ' ds ', s, 1.0d0/qs
dmax = max(dmax, 1.0d0/qs)
d = d + 1.0d0/qs
end do
@ -298,9 +299,9 @@ c
call ga_dgop(2, d , 1, '+')
c
if (ga_nodeid() .eq. 0) then
write(6,'(10x, i5, 2f15.10, 1p,d12.2,d12.2)')
write(luout,'(10x, i5, 2f15.10, 1p,d12.2,d12.2)')
& iter, dmax, d/dble(nloc), gamma_max
call util_flush(6)
call util_flush(luout)
end if
call ga_sync
c
@ -318,7 +319,7 @@ c
if (tt .le. nloc) then
s = iloc(ss)
t = iloc(tt)
* write(6,*) nloc, neven, set, pair, ss, tt, s, t
* write(luout,*) nloc, neven, set, pair, ss, tt, s, t
call ga_get(g_c, 1, nbf, s, s, c(1,1), 1)
call ga_get(g_sc, 1, nbf, s, s,sc(1,1), 1)
call ga_get(g_c, 1, nbf, t, t, c(1,2), 1)
@ -408,6 +409,7 @@ c
#include "geom.fh"
#include "bas.fh"
#include "util.fh"
#include "stdio.fh"
c
c Localize the nloc orbitals in iloc(*) by mixing with each other
c
@ -441,17 +443,17 @@ c
u21 = util_random(12345)
c
if (ga_nodeid() .eq. 0) then
write(6,9001)
write(luout,9001)
9001 format(/10x,' iter Max. dipole2 Mean dipole2 Converge'/
& 10x,' ---- ------------ ------------ --------')
call util_flush(6)
call util_flush(luout)
end if
c
dprev = 0.0d0
gamma_max = 0.0d0
call ga_distribution(g_uc(1), ga_nodeid(),
I ilo, ihi, jlo, jhi)
do iter = 1, 100
do iter = 1, 900
call ga_sync
nrot = 0
c
@ -496,9 +498,9 @@ c
call ga_dgop(2, d , 1, '+')
c
if (ga_nodeid() .eq. 0) then
write(6,'(10x, i5, 2f17.8, 1p,2d12.2)')
write(luout,'(10x, i5, 2f17.8, 1p,2d12.2)')
& iter, dmax, d/dble(nloc), gamma_max
call util_flush(6)
call util_flush(luout)
end if
call ga_sync
c
@ -521,7 +523,7 @@ cold do pair = 1+ga_nodeid(), neven/2, ga_nnodes()
c if g_uxc.s belongs to ga_nodeid, do
if(s.ge.jlo.and.s.le.jhi) then
t = iloc(tt)
* write(6,*) nloc, neven, set, pair, ss, tt, s, t
* write(luout,*) nloc, neven, set, pair, ss, tt, s, t
call ga_get(g_c,ilo, ihi, s, s, c(1,1), 1)
call ga_get(g_c,ilo, ihi, t, t, c(1,2), 1)
do x = 1, 4
@ -547,7 +549,7 @@ c
ast = ast + u12*u12 - 0.25d0*(u1-u2)**2
bst = bst + u12*(u1 - u2)
if (abs(u12-u21)/max(1.0d0,abs(u12)).gt.1d-8) then
write(6,*) ' U12, U21 ', u12, u21
write(luout,*) ' U12, U21 ', u12, u21
call errquit('bad u12', 0, UNKNOWN_ERR)
endif
end do
@ -584,13 +586,16 @@ c
call ga_sync
end do ! do set
end do ! do iter
if(ga_nodeid().eq.0) write(luout,*)
W ' WARNING: Foster-Boys solver not converged'
c
1000 continue
c
c Analyze localization of each mo
c
if (ga_nodeid() .eq. 0) then
write(6,*)
write(luout,*)
write(luout,*) ' Orbital Localization Analysis by atom'
do ss = 1, nloc
s = iloc(ss)
call ga_get(g_c, 1, nbf, s, s, c(1,1), 1)
@ -621,10 +626,10 @@ c
end if
end do
end do
write(6,9002) s, (list(a), pop(a), a=1,nlist)
write(luout,9002) s, (list(a), pop(a), a=1,nlist)
9002 format(i5, 100(2x,i4,'(',f5.2,')'))
end do
call util_flush(6)
call util_flush(luout)
end if
c
call ga_sync

1
src/property/GNUmakefile
View file

@ -65,6 +65,7 @@
localization_driver.o \
ibo_localization.o \
pm_localization.o \
fb_localization.o \
hnd_vec_write.o \
giao_b1_movecs_tools.o \
aor_r1_beta_anl_tools.o\

187
src/property/fb_localization.F Normal file
View file

@ -0,0 +1,187 @@
subroutine fb_localization(rtdb,geom,basis,
N nbf,nmo,nocc,nvir,
L ltyp,
G g_movecs)
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "msgids.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "bas.fh"
#include "util.fh"
#include "stdio.fh"
integer rtdb,geom,basis,
N nbf,nmo,nocc,nvir
character*3 ltyp
integer g_movecs
c
integer g_uc(4),g_smat,g_sc
integer nloc
integer l_sc, k_sc, l_c, k_c
integer l_iloc, k_iloc
integer l_dip(3), k_dip(3)
integer i,i_in
character*(15) pname
double precision sumd2,dx,dy,dz,rtmp
integer ga_create_atom_blocked
external ga_create_atom_blocked
logical debug
debug=.false.
c
c =================
c Boys localization
c =================
pname = 'fb_localization'
do i = 1,4
if (.not. ga_create(MT_DBL, nbf, nbf, 'uc',
$ nbf, 0, g_uc(i))) call errquit(pname//': uc'
& ,i,GA_ERR)
call ga_zero(g_uc(i))
end do
c dipole moment AO matrices -> uc(1) - uc(3)
call int_dip_ga(basis, basis, g_uc(1), g_uc(2), g_uc(3))
c calculate dipole-AO times C (mo-coeffs), store in uc(i)
if (.not. ga_create(MT_DBL, nbf, nmo, 'sc',
$ nbf, 0, g_sc)) call errquit(pname//': sc',0, GA_ERR)
do i = 1, 3
call ga_dgemm('n', 'n', nbf, nmo, nbf,
$ 1.0d0, g_uc(i), g_movecs, 0.0d0, g_sc)
call ga_copy_patch('n',g_sc,1,nbf,1,nmo,g_uc(i),1,nbf,1,nmo)
end do
c AO Overlap Matrix S:
g_smat = ga_create_atom_blocked(geom, basis, 'loc:smat')
call ga_zero(g_smat)
call int_1e_ga(basis, basis, g_smat, 'overlap', .false.)
c store S C in array uc(4)
call ga_dgemm('n', 'n', nbf, nmo, nbf,
$ 1.0d0, g_smat, g_movecs, 0.0d0, g_uc(4))
c
if (.not. ma_push_get(mt_dbl, 8*nbf, 'sc', l_sc, k_sc))
$ call errquit(pname//': ma for sc', 0, MA_ERR)
if (.not. ma_push_get(mt_dbl, 2*nbf, 'c', l_c, k_c))
$ call errquit(pname//': ma for c', 0, MA_ERR)
if (.not. ma_push_get(mt_int, nmo, 'iloc', l_iloc,
& k_iloc)) call errquit(pname//': loc:iloc', 1, MA_ERR)
if (debug) write(luout,*) 'MA c, sc complete'
c localize core and occupied orbitals:
if (ltyp.eq.'occ') then
do i = 1, nocc
int_mb(k_iloc - 1 + i) = i
c iloc(i) = i
enddo
nloc = nocc
else if (ltyp.eq.'vir') then
do i = nocc+1, nmo
int_mb(k_iloc - 1 - nocc + i) = i
c iloc(i-nocc) = i
end do
nloc = nmo - nocc
endif
c calculate orbital centroids and determine B2 = the sum of
c square distances of the centroids from the coordinate origin.
c The Boys localization maximizes B2.
c we can use g_sc as a temp array as it is not further used
if (.not. ga_destroy(g_sc)) call errquit(
& pname//': error destroying g_sc',0, GA_ERR)
if (.not. ga_create(MT_DBL, nmo, nmo, 'sc',
$ nbf, 0, g_sc)) call errquit(pname//': sc',0, GA_ERR)
do i = 1,3
if (.not. ma_push_get(mt_dbl, nmo, 'sc', l_dip(i), k_dip(i)))
$ call errquit(pname//': ma for dip', i, MA_ERR)
end do
do i = 1,3
call ga_zero(g_sc)
call ga_dgemm('t','n',nmo,nmo,nbf, 1.0d0, g_movecs,
& g_uc(i), 0.0d0, g_sc)
call ga_get_diagonal(g_sc, dbl_mb(k_dip(i)))
end do
do i = 3,1,-1
if (.not. ma_pop_stack(l_dip(i))) call errquit(
& pname//': error MA pop dip',i, MA_ERR)
end do
c jochen: comment:
c the Boys localization routine was already available
c in nwchem
call localizeFB(basis, dbl_mb(k_c), dbl_mb(k_sc),
$ nloc, int_mb(k_iloc), nbf, nmo, g_movecs, g_uc)
c calculate orbital centroids again and print information
do i = 1,3
if (.not. ma_push_get(mt_dbl, nmo, 'sc', l_dip(i), k_dip(i)))
$ call errquit(pname//': ma for dip', i, MA_ERR)
end do
do i = 1,3
call ga_zero(g_sc)
call ga_dgemm('t','n',nmo,nmo,nbf, 1.0d0, g_movecs,
& g_uc(i), 0.0d0, g_sc)
call ga_get_diagonal(g_sc, dbl_mb(k_dip(i)))
end do
if (ga_nodeid().eq.0) then
write (luout,'(/t6,a,t15,a,t27,a,t39,a,t49,a)')
& '#','<x>','<y>','<z>','<r>**2'
write (luout,'(1x,53(''-''))')
end if
sumd2 = 0.0d0
do i_in = 1,nloc
i=int_mb(k_iloc + i_in -1)
dx = dbl_mb(k_dip(1) + i-1)
dy = dbl_mb(k_dip(2) + i-1)
dz = dbl_mb(k_dip(3) + i-1)
rtmp = dx**2 + dy**2 + dz**2
sumd2 = sumd2 + rtmp
if (ga_nodeid().eq.0) then
write (luout,'(1x,i5,4f12.4)') i,dx,dy,dz,rtmp
end if
end do
c if (ga_nodeid().eq.0) then
c write (luout,'(1x,a,2f12.4)') 'B2 after localization, av:',
c & sumd2, sumd2/float(nloc)
c end if
if (.not. ma_chop_stack(l_dip(1))) call errquit(
& pname//': error MA pop dip',i, MA_ERR)
c clean up temp arrays:
if (.not. ga_destroy(g_sc)) call errquit(
& pname//': error destroying g_sc',0, GA_ERR)
if (.not. ma_chop_stack(l_sc)) call errquit(
& pname//': error MA pop sc',0, MA_ERR)
do i = 1, 4
if (.not. ga_destroy(g_uc(i)))
& call errquit(pname//': error destroying uc',i, GA_ERR)
end do
c smat not needed anymore
if (.not. ga_destroy(g_smat)) call errquit(
& pname//': error destroying g_smat',0, GA_ERR)
return
end

197
src/property/localization_driver.F
View file

@ -395,187 +395,36 @@ c =================
c Boys localization
c =================
if (nspin.gt.1) call
& errquit(pname//': Boys loc. not supported for open shell',
& nspin, INPUT_ERR)
do ispin = 1,nspin
ispin = 1
if (loc_opt.gt.0) call
& errquit(pname//': Boys loc. not supported for virtuals',
& loc_opt, INPUT_ERR)
do i = 1,4
if (.not. ga_create(MT_DBL, nbf, nbf, 'uc',
$ nbf, 0, g_uc(i))) call errquit(pname//': uc'
& ,i,GA_ERR)
call ga_zero(g_uc(i))
end do
if (debug) write(luout,*) 'g_uc created'
c dipole moment AO matrices -> uc(1) - uc(3)
call int_dip_ga(basis, basis, g_uc(1), g_uc(2), g_uc(3))
if (debug) write(luout,*) 'dipole done'
c calculate dipole-AO times C (mo-coeffs), store in uc(i)
if (.not. ga_create(MT_DBL, nbf, nmo, 'sc',
$ nbf, 0, g_sc)) call errquit(pname//': sc',0, GA_ERR)
do i = 1, 3
call ga_dgemm('n', 'n', nbf, nmo, nbf,
$ 1.0d0, g_uc(i), g_movecs(ispin), 0.0d0, g_sc)
call ga_copy_patch('n',g_sc,1,nbf,1,nmo,g_uc(i),1,nbf,1,nmo)
end do
if (debug) write(luout,*) 'g_uc complete'
c AO Overlap Matrix S:
g_smat = ga_create_atom_blocked(geom, basis, 'loc:smat')
call ga_zero(g_smat)
call int_1e_ga(basis, basis, g_smat, 'overlap', .false.)
if (debug) write(luout,*) 'smat done'
c store S C in array uc(4)
call ga_dgemm('n', 'n', nbf, nmo, nbf,
$ 1.0d0, g_smat, g_movecs(ispin), 0.0d0, g_uc(4))
if (debug) write(luout,*) 'g_uc(4) done'
c
if (.not. ma_push_get(mt_dbl, 8*nbf, 'sc', l_sc, k_sc))
$ call errquit(pname//': ma for sc', 0, MA_ERR)
if (.not. ma_push_get(mt_dbl, 8*nbf, 'c', l_c, k_c))
$ call errquit(pname//': ma for c', 0, MA_ERR)
if (.not. ma_push_get(mt_int, nmo, 'iloc', l_iloc,
& k_iloc)) call errquit(pname//': loc:iloc', 1, MA_ERR)
if (debug) write(luout,*) 'MA c, sc complete'
c localize core and occupied orbitals:
do i = 1, nclosed(1)
int_mb(k_iloc - 1 + i) = i
c iloc(i) = i
end do
nloc = nclosed(1)
c calculate orbital centroids and determine B2 = the sum of
c square distances of the centroids from the coordinate origin.
c The Boys localization maximizes B2.
c we can use g_sc as a temp array as it is not further used
if (.not. ga_destroy(g_sc)) call errquit(
& pname//': error destroying g_sc',0, GA_ERR)
if (.not. ga_create(MT_DBL, nmo, nmo, 'sc',
$ nbf, 0, g_sc)) call errquit(pname//': sc',0, GA_ERR)
do i = 1,3
if (.not. ma_push_get(mt_dbl, nmo, 'sc', l_dip(i), k_dip(i)))
$ call errquit(pname//': ma for dip', i, MA_ERR)
end do
do i = 1,3
call ga_zero(g_sc)
call ga_dgemm('t','n',nmo,nmo,nbf, 1.0d0, g_movecs(ispin),
& g_uc(i), 0.0d0, g_sc)
call ga_get_diagonal(g_sc, dbl_mb(k_dip(i)))
end do
if (ga_nodeid().eq.0) then
write (luout,'(/t6,a,t15,a,t27,a,t39,a,t49,a)')
& '#','<x>','<y>','<z>','<r>**2'
write (luout,'(1x,53(''-''))')
end if
sumd2 = 0.0d0
do i = 1,nmo
dx = dbl_mb(k_dip(1) + i-1)
dy = dbl_mb(k_dip(2) + i-1)
dz = dbl_mb(k_dip(3) + i-1)
rtmp = dx**2 + dy**2 + dz**2
sumd2 = sumd2 + rtmp
if (ga_nodeid().eq.0) then
write (luout,'(1x,i5,4f12.4)') i,dx,dy,dz,rtmp
if (nspin.eq.1) then
nocc = nclosed(ispin)
else if (nspin.eq.2) then
nocc = nopen(ispin)
end if
end do
c if (ga_nodeid().eq.0) then
c write (luout,'(1x,a,2f12.4)') 'B2 before localization, av:',
c & sumd2, sumd2/float(nloc)
c end if
nvir = nvirt(ispin)
do i = 3,1,-1
if (.not. ma_pop_stack(l_dip(i))) call errquit(
& pname//': error MA pop dip',i, MA_ERR)
end do
c jochen: comment:
c the Boys localization routine was already available
c in nwchem
call localizeFB(basis, dbl_mb(k_c), dbl_mb(k_sc),
$ nloc, int_mb(k_iloc), nbf, nmo, g_movecs(ispin), g_uc)
c calculate orbital centroids again and print information
do i = 1,3
if (.not. ma_push_get(mt_dbl, nmo, 'sc', l_dip(i), k_dip(i)))
$ call errquit(pname//': ma for dip', i, MA_ERR)
end do
do i = 1,3
call ga_zero(g_sc)
call ga_dgemm('t','n',nmo,nmo,nbf, 1.0d0, g_movecs(ispin),
& g_uc(i), 0.0d0, g_sc)
call ga_get_diagonal(g_sc, dbl_mb(k_dip(i)))
end do
if (ga_nodeid().eq.0) then
write (luout,'(/t6,a,t15,a,t27,a,t39,a,t49,a)')
& '#','<x>','<y>','<z>','<r>**2'
write (luout,'(1x,53(''-''))')
end if
sumd2 = 0.0d0
do i = 1,nmo
dx = dbl_mb(k_dip(1) + i-1)
dy = dbl_mb(k_dip(2) + i-1)
dz = dbl_mb(k_dip(3) + i-1)
rtmp = dx**2 + dy**2 + dz**2
sumd2 = sumd2 + rtmp
if (ga_nodeid().eq.0) then
write (luout,'(1x,i5,4f12.4)') i,dx,dy,dz,rtmp
if (nspin.gt.1 .and. master) then
write(luout,"(/1x,6('=')/1x,'Spin ',i1/1x,6('='))") ispin
end if
end do
c if (ga_nodeid().eq.0) then
c write (luout,'(1x,a,2f12.4)') 'B2 after localization, av:',
c & sumd2, sumd2/float(nloc)
c end if
do i = 3,1,-1
if (.not. ma_pop_stack(l_dip(i))) call errquit(
& pname//': error MA pop dip',i, MA_ERR)
end do
if (loc_opt.eq.0 .or. loc_opt.eq.2) then
c clean up temp arrays:
if (.not. ga_destroy(g_sc)) call errquit(
& pname//': error destroying g_sc',0, GA_ERR)
if (.not. ma_pop_stack(l_c)) call errquit(
& pname//': error MA pop c',0, MA_ERR)
if (.not. ma_pop_stack(l_sc)) call errquit(
& pname//': error MA pop sc',0, MA_ERR)
do i = 1, 4
if (.not. ga_destroy(g_uc(i)))
& call errquit(pname//': error destroying uc',i, GA_ERR)
end do
c smat not needed anymore
if (.not. ga_destroy(g_smat)) call errquit(
& pname//': error destroying g_smat',0, GA_ERR)
ltyp = 'occ'
call fb_localization(rtdb,geom,basis,
N nbf,nmo,nocc,nvir,
L ltyp,
G g_movecs(ispin))
end if
if (loc_opt.eq.1 .or. loc_opt.eq.2) then
ltyp = 'vir'
call fb_localization(rtdb,geom,basis,
N nbf,nmo,nocc,nvir,
L ltyp,
G g_movecs(ispin))
end if
enddo
c ===============================================================