From 0a1903a063faa5b873f250f006d4a4a1dae7e11e Mon Sep 17 00:00:00 2001 From: edoapra Date: Wed, 25 Mar 2026 16:20:14 -0700 Subject: [PATCH 1/2] increased no. of iterations for Foster-Boys --- src/ddscf/localize.F | 37 +++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/src/ddscf/localize.F b/src/ddscf/localize.F index 622026506f..e38da3de6c 100644 --- a/src/ddscf/localize.F +++ b/src/ddscf/localize.F @@ -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 From 2e59c51500260e8523bbac14a019b0ff43a231c3 Mon Sep 17 00:00:00 2001 From: edoapra Date: Wed, 25 Mar 2026 16:21:41 -0700 Subject: [PATCH 2/2] extracted Foster-Boys interface to allow localization of virtual orbitals --- src/property/GNUmakefile | 1 + src/property/fb_localization.F | 187 +++++++++++++++++++++++++++ src/property/localization_driver.F | 197 ++++------------------------- 3 files changed, 211 insertions(+), 174 deletions(-) create mode 100644 src/property/fb_localization.F diff --git a/src/property/GNUmakefile b/src/property/GNUmakefile index 52145a27dc..be24fcbcac 100644 --- a/src/property/GNUmakefile +++ b/src/property/GNUmakefile @@ -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\ diff --git a/src/property/fb_localization.F b/src/property/fb_localization.F new file mode 100644 index 0000000000..6171bc87c5 --- /dev/null +++ b/src/property/fb_localization.F @@ -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)') + & '#','','','','**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 diff --git a/src/property/localization_driver.F b/src/property/localization_driver.F index b8b11ef6d5..f7d83a0920 100644 --- a/src/property/localization_driver.F +++ b/src/property/localization_driver.F @@ -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)') - & '#','','','','**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)') - & '#','','','','**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 ===============================================================