diff --git a/QA/nwparse.pl b/QA/nwparse.pl index 73cb0c66e0..974bafd720 100644 --- a/QA/nwparse.pl +++ b/QA/nwparse.pl @@ -402,7 +402,7 @@ foreach $filename (@FILES_TO_PARSE) { } } - if (/Excitation energy/ || /Rotatory /) { + if (/Excitation energy/ || /Rotatory / || /IBO loc: largest element in /) { if ($debug) {print "\ndebug: $_";} @line_tokens = split(' '); $num_line_tokens = @line_tokens; diff --git a/QA/tests/localize-ibo-aa/localize-ibo-aa.nw b/QA/tests/localize-ibo-aa/localize-ibo-aa.nw index 049887f2a7..3a1b422a5c 100644 --- a/QA/tests/localize-ibo-aa/localize-ibo-aa.nw +++ b/QA/tests/localize-ibo-aa/localize-ibo-aa.nw @@ -71,7 +71,7 @@ convergence energy 1e-10 end property - localization ibo + localization ibo 2 end task scf property diff --git a/QA/tests/localize-ibo-aa/localize-ibo-aa.out b/QA/tests/localize-ibo-aa/localize-ibo-aa.out index ee86c74cb0..adfec65d11 100644 --- a/QA/tests/localize-ibo-aa/localize-ibo-aa.out +++ b/QA/tests/localize-ibo-aa/localize-ibo-aa.out @@ -1,4 +1,4 @@ -running on 4 processors +running on 6 processors argument 1 = acrylic-acid.nw @@ -77,12 +77,14 @@ convergence energy 1e-10 end property - localization ibo + localization ibo 2 end +#task scf energy task scf property + ================================================================================ @@ -90,7 +92,7 @@ task scf property - Northwest Computational Chemistry Package (NWChem) 7.0.1 + Northwest Computational Chemistry Package (NWChem) 7.2.0 -------------------------------------------------------- @@ -98,7 +100,7 @@ task scf property Pacific Northwest National Laboratory Richland, WA 99352 - Copyright (c) 1994-2020 + Copyright (c) 1994-2022 Pacific Northwest National Laboratory Battelle Memorial Institute @@ -123,21 +125,21 @@ task scf property Job information --------------- - hostname = ja31 - program = /home/workspace/jochena/nwchem/github-fork-ibo-devel/bin/LINUX64/nwchem - date = Thu May 5 15:37:05 2022 + hostname = ja04 + program = /home/workspace/jochena/nwchem/github-fork-ibovir-devel/bin/LINUX64/nwchem + date = Tue Feb 14 12:19:49 2023 - compiled = Thu_May_05_15:36:59_2022 - source = /home/workspace/jochena/nwchem/github-fork-ibo-devel - nwchem branch = 7.0.0 - nwchem revision = nwchem_on_git-3953-g107e8544ce - ga revision = 5.8.1 - use scalapack = T + compiled = Tue_Feb_14_12:19:44_2023 + source = /home/workspace/jochena/nwchem/github-fork-ibovir-devel + nwchem branch = 7.2.0 + nwchem revision = nwchem_on_git-4622-g8f291f7cfd + ga revision = 5.8.0 + use scalapack = F input = acrylic-acid.nw prefix = testjob. data base = ./testjob.db status = startup - nproc = 3 + nproc = 5 time left = -1s @@ -145,10 +147,10 @@ task scf property Memory information ------------------ - heap = 6553596 doubles = 50.0 Mbytes - stack = 6553601 doubles = 50.0 Mbytes + heap = 6553598 doubles = 50.0 Mbytes + stack = 6553595 doubles = 50.0 Mbytes global = 13107200 doubles = 100.0 Mbytes (distinct from heap & stack) - total = 26214397 doubles = 200.0 Mbytes + total = 26214393 doubles = 200.0 Mbytes verify = yes hardfail = no @@ -265,7 +267,7 @@ task scf property library name resolved from: environment library file name is: < - /home/workspace/jochena/nwchem/github-fork-ibo-devel/src/basis/libraries/> + /home/workspace/jochena/nwchem/github-fork-ibovir-devel/src/basis/libraries/> @@ -278,7 +280,7 @@ task scf property library name resolved from: environment library file name is: < - /home/workspace/jochena/nwchem/github-fork-ibo-devel/src/basis/libraries/> + /home/workspace/jochena/nwchem/github-fork-ibovir-devel/src/basis/libraries/> Basis "iao basis" -> "" (spherical) ----- @@ -502,22 +504,22 @@ task scf property iter energy gnorm gmax time ----- ------------------- --------- --------- -------- - 1 -265.3310324443 1.42D+00 2.92D-01 0.2 - 2 -265.4360577045 4.13D-01 8.78D-02 0.4 - 3 -265.4514376314 3.64D-02 9.22D-03 0.8 - 4 -265.4516512505 1.03D-03 2.34D-04 1.3 - 5 -265.4516513710 7.80D-06 1.92D-06 1.9 + 1 -265.3310324443 1.42D+00 2.92D-01 0.1 + 2 -265.4360577045 4.13D-01 8.78D-02 0.2 + 3 -265.4514376314 3.64D-02 9.22D-03 0.4 + 4 -265.4516512505 1.03D-03 2.34D-04 0.6 + 5 -265.4516513711 7.80D-06 1.92D-06 0.9 Final RHF results ------------------ - Total SCF energy = -265.451651371048 - One-electron energy = -682.442839763839 + Total SCF energy = -265.451651371050 + One-electron energy = -682.442839763841 Two-electron energy = 255.473969802116 Nuclear repulsion energy = 161.517218590675 - Time for solution = 1.8s + Time for solution = 0.9s Final eigenvalues @@ -558,7 +560,7 @@ task scf property ------------------------------------- Vector 6 Occ=2.000000D+00 E=-1.441125D+00 - MO Center= -1.1D+00, 2.6D-02, -3.4D-17, r^2= 1.2D+00 + MO Center= -1.1D+00, 2.6D-02, -2.1D-17, r^2= 1.2D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 58 0.413566 5 O s 30 0.293845 3 C s @@ -566,7 +568,7 @@ task scf property 57 0.208468 5 O s 45 0.150994 4 O s Vector 7 Occ=2.000000D+00 E=-1.338769D+00 - MO Center= -9.1D-01, -4.1D-01, -2.6D-17, r^2= 1.4D+00 + MO Center= -9.1D-01, -4.1D-01, 6.6D-17, r^2= 1.4D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 44 0.424994 4 O s 58 -0.352816 5 O s @@ -575,7 +577,7 @@ task scf property 57 -0.176450 5 O s Vector 8 Occ=2.000000D+00 E=-1.070923D+00 - MO Center= 1.1D+00, 2.8D-01, -4.7D-17, r^2= 1.5D+00 + MO Center= 1.1D+00, 2.8D-01, 3.4D-17, r^2= 1.5D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 16 0.400200 2 C s 2 0.337172 1 C s @@ -583,7 +585,7 @@ task scf property 17 0.160365 2 C s Vector 9 Occ=2.000000D+00 E=-8.974613D-01 - MO Center= 2.3D-01, 2.6D-01, -9.8D-17, r^2= 3.3D+00 + MO Center= 2.3D-01, 2.6D-01, -1.0D-16, r^2= 3.3D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 30 -0.305983 3 C s 2 0.303431 1 C s @@ -592,7 +594,7 @@ task scf property 1 -0.152955 1 C s Vector 10 Occ=2.000000D+00 E=-7.642057D-01 - MO Center= 1.6D-01, 3.0D-01, 1.4D-16, r^2= 4.1D+00 + MO Center= 1.6D-01, 3.0D-01, 5.3D-16, r^2= 4.1D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 60 0.316925 5 O px 16 -0.241192 2 C s @@ -602,7 +604,7 @@ task scf property 4 0.156757 1 C px 71 0.155748 6 H s Vector 11 Occ=2.000000D+00 E=-7.114488D-01 - MO Center= 2.7D-01, 9.8D-03, -1.4D-16, r^2= 3.1D+00 + MO Center= 2.7D-01, 9.8D-03, -3.3D-16, r^2= 3.1D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 19 0.248395 2 C py 81 0.212079 8 H s @@ -612,7 +614,7 @@ task scf property 44 0.168207 4 O s 46 -0.165171 4 O px Vector 12 Occ=2.000000D+00 E=-6.942998D-01 - MO Center= -5.7D-01, -1.8D-01, 7.4D-17, r^2= 3.0D+00 + MO Center= -5.7D-01, -1.8D-01, -2.1D-16, r^2= 3.0D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 61 0.332130 5 O py 47 -0.285481 4 O py @@ -622,7 +624,7 @@ task scf property 33 0.150516 3 C py Vector 13 Occ=2.000000D+00 E=-6.290129D-01 - MO Center= -9.9D-01, 1.2D-01, 2.6D-16, r^2= 1.6D+00 + MO Center= -9.9D-01, 1.2D-01, -3.2D-16, r^2= 1.6D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 62 0.438625 5 O pz 65 0.308876 5 O pz @@ -630,7 +632,7 @@ task scf property 51 0.167682 4 O pz 37 0.158757 3 C pz Vector 14 Occ=2.000000D+00 E=-6.234118D-01 - MO Center= 8.0D-01, 1.3D-02, -1.6D-16, r^2= 3.8D+00 + MO Center= 8.0D-01, 1.3D-02, 4.5D-16, r^2= 3.8D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 4 0.334286 1 C px 76 0.262842 7 H s @@ -639,7 +641,7 @@ task scf property 18 -0.161739 2 C px Vector 15 Occ=2.000000D+00 E=-5.982438D-01 - MO Center= -2.3D-02, 2.4D-01, -5.0D-17, r^2= 3.5D+00 + MO Center= -2.3D-02, 2.4D-01, -3.8D-17, r^2= 3.5D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 61 0.305744 5 O py 18 0.273635 2 C px @@ -649,7 +651,7 @@ task scf property 33 -0.180575 3 C py 71 -0.170295 6 H s Vector 16 Occ=2.000000D+00 E=-5.496397D-01 - MO Center= 1.1D+00, 6.0D-02, 2.5D-17, r^2= 3.3D+00 + MO Center= 1.1D+00, 6.0D-02, 3.0D-16, r^2= 3.3D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 5 0.298431 1 C py 71 -0.271981 6 H s @@ -658,7 +660,7 @@ task scf property 33 0.163556 3 C py 61 -0.150120 5 O py Vector 17 Occ=2.000000D+00 E=-4.889073D-01 - MO Center= -9.0D-01, -2.5D-01, -6.7D-17, r^2= 2.0D+00 + MO Center= -9.0D-01, -2.5D-01, -5.3D-16, r^2= 2.0D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 62 -0.430815 5 O pz 48 0.407766 4 O pz @@ -666,7 +668,7 @@ task scf property 34 0.166236 3 C pz Vector 18 Occ=2.000000D+00 E=-4.582158D-01 - MO Center= -5.2D-01, -7.3D-01, -9.6D-17, r^2= 1.9D+00 + MO Center= -5.2D-01, -7.3D-01, -1.6D-16, r^2= 1.9D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 46 0.549426 4 O px 49 0.417524 4 O px @@ -674,7 +676,7 @@ task scf property 64 -0.174821 5 O py 18 0.151160 2 C px Vector 19 Occ=2.000000D+00 E=-4.009704D-01 - MO Center= 1.1D+00, 1.5D-01, -2.0D-16, r^2= 2.0D+00 + MO Center= 1.1D+00, 1.5D-01, -4.1D-16, r^2= 2.0D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 20 0.384445 2 C pz 6 0.352495 1 C pz @@ -682,7 +684,7 @@ task scf property 48 -0.214439 4 O pz 51 -0.173233 4 O pz Vector 20 Occ=0.000000D+00 E= 9.067747D-02 - MO Center= 8.9D-01, 1.5D-02, 2.7D-17, r^2= 3.0D+00 + MO Center= 8.9D-01, 1.5D-02, 9.9D-17, r^2= 3.0D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 9 0.621145 1 C pz 23 -0.459696 2 C pz @@ -692,7 +694,7 @@ task scf property 65 0.155474 5 O pz Vector 21 Occ=0.000000D+00 E= 1.933962D-01 - MO Center= -2.6D+00, 1.3D-01, -3.0D-16, r^2= 2.6D+00 + MO Center= -2.6D+00, 1.3D-01, 2.0D-16, r^2= 2.6D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 87 1.532433 9 H s 17 -0.762684 2 C s @@ -702,7 +704,7 @@ task scf property 60 0.181785 5 O px 3 0.180050 1 C s Vector 22 Occ=0.000000D+00 E= 1.990748D-01 - MO Center= 2.4D+00, 1.2D+00, -1.1D-16, r^2= 4.4D+00 + MO Center= 2.4D+00, 1.2D+00, 5.2D-17, r^2= 4.4D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 3 1.857859 1 C s 77 -1.791874 7 H s @@ -712,7 +714,7 @@ task scf property 8 0.218519 1 C py 35 0.169789 3 C px Vector 23 Occ=0.000000D+00 E= 2.314259D-01 - MO Center= 1.7D+00, -2.0D-01, -9.3D-17, r^2= 6.0D+00 + MO Center= 1.7D+00, -2.0D-01, 9.0D-16, r^2= 6.0D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 3 1.936960 1 C s 72 -1.805329 6 H s @@ -721,7 +723,7 @@ task scf property 8 -0.489914 1 C py 5 -0.154259 1 C py Vector 24 Occ=0.000000D+00 E= 2.626515D-01 - MO Center= 1.3D-01, -2.5D-02, -4.7D-16, r^2= 2.9D+00 + MO Center= 1.3D-01, -2.5D-02, 7.6D-16, r^2= 2.9D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 37 -0.804941 3 C pz 23 0.756689 2 C pz @@ -731,7 +733,7 @@ task scf property 62 0.158924 5 O pz Vector 25 Occ=0.000000D+00 E= 2.688800D-01 - MO Center= 1.9D+00, 3.8D-01, 5.4D-16, r^2= 6.9D+00 + MO Center= 1.9D+00, 3.8D-01, -1.8D-15, r^2= 6.9D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 77 -1.706163 7 H s 72 1.689090 6 H s @@ -741,7 +743,7 @@ task scf property 87 -0.319009 9 H s 3 -0.288918 1 C s Vector 26 Occ=0.000000D+00 E= 3.855495D-01 - MO Center= -4.6D-01, -2.3D-01, 1.1D-16, r^2= 3.7D+00 + MO Center= -4.6D-01, -2.3D-01, -1.7D-16, r^2= 3.7D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 31 4.026240 3 C s 45 -1.439431 4 O s @@ -751,7 +753,7 @@ task scf property 50 -0.560923 4 O py 82 -0.528079 8 H s Vector 27 Occ=0.000000D+00 E= 4.475195D-01 - MO Center= 5.8D-01, 2.5D-01, -9.4D-17, r^2= 3.9D+00 + MO Center= 5.8D-01, 2.5D-01, 3.4D-16, r^2= 3.9D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 21 3.533345 2 C px 3 -2.731556 1 C s @@ -761,7 +763,7 @@ task scf property 7 0.878572 1 C px 30 0.524028 3 C s Vector 28 Occ=0.000000D+00 E= 4.523376D-01 - MO Center= 1.5D+00, -1.4D-01, 9.3D-17, r^2= 4.3D+00 + MO Center= 1.5D+00, -1.4D-01, -3.7D-16, r^2= 4.3D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 17 6.048068 2 C s 3 -4.672145 1 C s @@ -771,7 +773,7 @@ task scf property 59 -0.983223 5 O s 36 0.975052 3 C py Vector 29 Occ=0.000000D+00 E= 4.789847D-01 - MO Center= 2.6D-01, 3.3D-01, -1.0D-16, r^2= 3.8D+00 + MO Center= 2.6D-01, 3.3D-01, 2.1D-18, r^2= 3.8D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 17 2.858248 2 C s 3 -2.090889 1 C s @@ -781,7 +783,7 @@ task scf property 35 0.814942 3 C px 63 0.614144 5 O px Vector 30 Occ=0.000000D+00 E= 5.702926D-01 - MO Center= 1.4D+00, 5.8D-01, 2.3D-15, r^2= 2.8D+00 + MO Center= 1.4D+00, 5.8D-01, -1.9D-15, r^2= 2.8D+00 Bfn. Coefficient Atom+Function Bfn. Coefficient Atom+Function ----- ------------ --------------- ----- ------------ --------------- 17 2.332617 2 C s 31 -1.303425 3 C s @@ -831,7 +833,7 @@ task scf property 2 1 1 0 0.077897 0.000000 0.000002 2 1 0 1 0.000000 0.000000 0.000000 2 0 2 0 -24.880694 0.000000 90.697964 - 2 0 1 1 0.000000 0.000000 0.000000 + 2 0 1 1 -0.000000 0.000000 0.000000 2 0 0 2 -22.450650 0.000000 0.000000 @@ -842,19 +844,20 @@ task scf property *** ORBITAL LOCALIZATION *** **************************** - IAO-IBO localization - -------------------- + IAO-IBO localization + -------------------- + (occ and virt orbitals) - IBO loc: largest element in C(iao, trans) S C(iao) -1: 0.000000000000 + IBO loc: largest element in C(iao,T) S C(iao) -1: 0.00000000 Significant deviations from zero may indicate elevated numerical noise in the IAO generation - IBO loc: largest element of C(MO, trans) C(MO) -1: 0.000000000000 + IBO loc: largest element in C(MO,T) C(MO) -1: 0.00000000 should be zero, for CMOs in the IAO basis iter Max. delocal Mean delocal Converge - ---- ------------ ------------ --------- + ---- ------------ ------------ -------- 1 6.4744549430 3.2718681346 0.00D+00 2 2.9502374379 1.7407413313 7.84D-01 3 2.1535052593 1.5156908043 6.80D-01 @@ -863,10 +866,10 @@ task scf property 6 2.1555070209 1.5083589529 8.63D-05 7 2.1555070209 1.5083589538 6.36D-06 8 2.1555070209 1.5083589538 9.22D-08 - 9 2.1555070209 1.5083589538 0.00D+00 + 9 2.1555070209 1.5083589538 3.73D-09 - IBO loc: largest element of C(MO, trans) C(MO) -1: 0.000000000000 - should be zero, for LMOs in the IAO basis + IBO loc: largest element in C(MO,T) C(MO) -1: 0.00000000 + should be zero, for IBOs in the IAO basis IAO-IBO localized orbitals @@ -890,9 +893,59 @@ task scf property 18 -0.507592 2.000 4( 0.93) 3( 0.04) 5( 0.01) 2( 0.01) 19 -0.424380 2.000 2( 0.53) 1( 0.43) 3( 0.04) + + IBO localization (occ): IBOs will be stored + in file locorb.movecs, number + 1 to 19 + + + IBO loc: largest element in C(iao,T) S C(iao) -1: 0.00000000 + Significant deviations from zero may indicate + elevated numerical noise in the IAO generation + + non-zero singular values: 10 + + IBO loc: largest element in C(MO,T) C(MO) -1: 0.00000000 + should be zero, for CMOs in the IAO basis + + iter Max. delocal Mean delocal Converge + ---- ------------ ------------ -------- + 1 5.7963239473 4.2732750386 0.00D+00 + 2 2.4941110119 2.1820228090 7.23D-01 + 3 2.1364906389 1.9920797446 4.32D-01 + 4 2.1348072188 1.9904017968 1.39D-01 + 5 2.1348013497 1.9903805756 9.02D-03 + 6 2.1348017491 1.9903804736 7.84D-04 + 7 2.1348017624 1.9903804751 5.13D-05 + 8 2.1348017624 1.9903804751 3.41D-06 + 9 2.1348017624 1.9903804751 2.26D-07 + 10 2.1348017624 1.9903804751 1.34D-08 + 11 2.1348017624 1.9903804751 5.27D-09 + + IBO loc: largest element in C(MO,T) C(MO) -1: 0.00000000 + should be zero, for IBOs in the IAO basis + + IAO-IBO localized orbitals + + 1 0.209314 0.000 1( 0.55) 2( 0.44) + 2 0.251712 0.000 3( 0.64) 4( 0.25) 5( 0.07) 2( 0.03) + 3 0.441767 0.000 9( 0.69) 5( 0.30) + 4 0.571110 0.000 7( 0.57) 1( 0.42) + 5 0.576430 0.000 6( 0.58) 1( 0.41) + 6 0.588168 0.000 8( 0.57) 2( 0.42) + 7 0.612417 0.000 3( 0.61) 5( 0.34) 4( 0.04) + 8 0.678788 0.000 2( 0.50) 3( 0.47) 4( 0.02) + 9 0.845663 0.000 3( 0.59) 4( 0.40) + 10 0.862369 0.000 1( 0.50) 2( 0.50) + + + IBO localization (vir): IBOs will be stored + in file locorb.movecs, number + 20 to 29 + Exiting Localization driver routine - Task times cpu: 2.0s wall: 2.1s + Task times cpu: 1.0s wall: 1.2s NWChem Input Module @@ -914,7 +967,7 @@ MA usage statistics: current number of blocks 0 0 maximum number of blocks 22 14 current total bytes 0 0 - maximum total bytes 139392 22511336 + maximum total bytes 139376 22511288 maximum total K-bytes 140 22512 maximum total M-bytes 1 23 @@ -954,10 +1007,10 @@ MA usage statistics: ------- E. Apra, E. J. Bylaska, N. Govind, K. Kowalski, M. Valiev, D. Mejia-Rodriguez, A. Kunitsa, N. P. Bauman, A. Panyala, W. A. de Jong, T. P. Straatsma, - H. J. J. van Dam, D. Wang, T. L. Windus, J. Hammond, J. Autschbach, + H. J. J. van Dam, D. Wang, T. L. Windus, J. Hammond, J. Autschbach, A. Woods, K. Bhaskaran-Nair, J. Brabec, K. Lopata, S. A. Fischer, S. Krishnamoorthy, M. Jacquelin, W. Ma, M. Klemm, O. Villa, Y. Chen, V. Anisimov, F. Aquino, - S. Hirata, M. T. Hackler, Eric Hermes, L. Jensen, J. E. Moore, J. C. Becca, + S. Hirata, M. T. Hackler, E. Hermes, L. Jensen, J. E. Moore, J. C. Becca, V. Konjkov, T. Risthaus, M. Malagoli, A. Marenich, A. Otero-de-la-Roza, J. Mullin, P. Nichols, R. Peverati, J. Pittner, Y. Zhao, P.-D. Fan, A. Fonari, M. J. Williamson, R. J. Harrison, J. R. Rehr, M. Dupuis, @@ -972,4 +1025,4 @@ MA usage statistics: K. Glaesemann, G. Sandrone, M. Stave, H. Taylor, G. Thomas, J. H. van Lenthe, A. T. Wong, Z. Zhang. - Total times cpu: 2.0s wall: 2.1s + Total times cpu: 1.0s wall: 1.2s diff --git a/src/ddscf/localize.F b/src/ddscf/localize.F index 6e42a46c1b..8b0862c79a 100644 --- a/src/ddscf/localize.F +++ b/src/ddscf/localize.F @@ -276,9 +276,9 @@ c gamma_tol = 1d-10 c if (ga_nodeid() .eq. 0) then - write(6,2) - 2 format(/10x,' iter Max. delocal Mean delocal Converge'/ - $ 10x,' ---- ------------ ------------ ---------') + write(6,9001) + 9001 format(/10x,' iter Max. delocal Mean delocal Converge'/ + & 10x,' ---- ------------ ------------ --------') call util_flush(6) end if c @@ -323,9 +323,9 @@ c call ga_dgop(2, d , 1, '+') c if (ga_nodeid() .eq. 0) then - write(6,1) iter, dmax, d/dble(nloc), gamma_max - 1 format(10x, i5, 2f15.10, 1p,d12.2,d12.2) - call util_flush(6) + write(6,'(10x, i5, 2f15.10, 1p,d12.2,d12.2)') + & iter, dmax, d/dble(nloc), gamma_max + call util_flush(6) end if call ga_sync c @@ -439,8 +439,8 @@ c end if end do end do - write(6,77) s, (list(a), pop(a), a=1,nlist) - 77 format(i5, 100(2x,i4,'(',f5.2,')')) + write(6,9002) s, (list(a), pop(a), a=1,nlist) + 9002 format(i5, 100(2x,i4,'(',f5.2,')')) end do call util_flush(6) end if @@ -498,10 +498,10 @@ c u21 = util_random(12345) c if (ga_nodeid() .eq. 0) then - write(6,2) - 2 format(/10x,' iter Max. dipole2 Mean dipole2 Converge'/ - $ 10x,' ---- ------------ ------------ ---------') - call util_flush(6) + write(6,9001) + 9001 format(/10x,' iter Max. dipole2 Mean dipole2 Converge'/ + & 10x,' ---- ------------ ------------ --------') + call util_flush(6) end if c dprev = 0.0d0 @@ -544,9 +544,9 @@ c call ga_dgop(2, d , 1, '+') c if (ga_nodeid() .eq. 0) then - write(6,1) iter, dmax, d/dble(nloc), gamma_max - 1 format(10x, i5, 2f17.8, 1p,2d12.2) - call util_flush(6) + write(6,'(10x, i5, 2f17.8, 1p,2d12.2)') + & iter, dmax, d/dble(nloc), gamma_max + call util_flush(6) end if call ga_sync c @@ -663,8 +663,8 @@ c end if end do end do - write(6,77) s, (list(a), pop(a), a=1,nlist) - 77 format(i5, 100(2x,i4,'(',f5.2,')')) + write(6,9002) s, (list(a), pop(a), a=1,nlist) + 9002 format(i5, 100(2x,i4,'(',f5.2,')')) end do call util_flush(6) end if @@ -674,16 +674,19 @@ c end c c ================================================================= -c - subroutine localizeIBO(basis, minbas, - & c, nloc, iloc, nbf, nmo, - & g_c, g_smat, g_iao, mnbf, eval, occ) c -c IAO-IBO localization -c - + subroutine localizeIBO(minbas, c, g_c, nloc, nbf, mnbf, + & natoms) +c +c IBO localization. Similar to Pipek-Mezey, but in an orthonormal +c IAO basis. Analysis of the resulting IBOs is done in the calling +c routine. +c +c written by Jochen Autschbach, jochena@buffalo.edu, +c adapted from routine localizePM + implicit none - + #include "errquit.fh" #include "nwc_const.fh" #include "mafdecls.fh" @@ -692,136 +695,92 @@ c #include "bas.fh" #include "util.fh" #include "stdio.fh" -c -c Localize the nloc orbitals in iloc(*) by mixing with each other -c - integer basis, minbas, nloc, iloc(*), nbf, nmo, mnbf - double precision c(mnbf, 2), eval(nbf), occ(nbf) - integer g_c, g_smat, g_iao - integer maxat, nlist -*............................ these should be dynamically allocated ? - parameter (maxat = nw_max_atom) - integer list(maxat) - double precision pop(maxat) - integer iter, ss, s, tt, t, a, u, bflo, bfhi, natoms, geom +c subroutine arguments: + + integer minbas ! minimal basis for IAOs (in) + integer nloc ! number of orbitals to localize (in) + integer nbf ! size of AO basis (in) + integer mnbf ! size of minimal basis (in) + integer natoms ! number of atoms (in) + double precision c(nbf, 2) ! memory used for MO transformations + integer g_c ! orbital coefficients in IAO basis (inout) + +c local variables: + + integer nlist + integer iter, s, t, a, u, bflo, bfhi, geom double precision ast, bst, qast, qat, qas, gamma, cosg, sing, d, & qs, dprev, tol, dmax, gamma_tol, gamma_max, tmp double precision minval, swap integer nrot, set, pair, neven - integer ntmp1, ntmp2 logical master, debug - integer imo, jmo, kmo - - integer g_tmp1, g_tmp2, g_cibocc - - character*(11) pname - - character*(256) lmotrans - logical file_write_ga - external file_write_ga + integer g_tmp + character*(12) pname c ----------------------------------------------------------------- - pname = 'iao-ibo loc' + pname = 'iao-ibo iter' master = (ga_nodeid() .eq. 0) debug = .false. ! true, during development if (debug.and.master) write(luout,*) 'entering '//pname -c - if (.not. bas_geom(basis, geom)) call errquit - $ (pname//': basis ', 0, BASIS_ERR) + +c sanity checks, so we don't crash later: + + if (.not. bas_geom(minbas, geom)) call errquit + & (pname//': minbas ', 0, BASIS_ERR) if (.not. geom_ncent(geom, natoms)) call errquit - $ (pname//': geom',0, GEOM_ERR) -c - if (natoms.gt.maxat) call errquit - & (pname//': maxat too small ', 911, UNKNOWN_ERR) + & (pname//': geom',0, GEOM_ERR) + if (mnbf>nbf) call errquit + & (pname//': mnbf>nbf', 0, UNKNOWN_ERR) + if (nloc>nbf) call errquit + & (pname//': nloc>nbf', 0, UNKNOWN_ERR) + if (nloc>mnbf) call errquit + & (pname//': nloc>mnbf', 0, UNKNOWN_ERR) -c note sure if we need this for sure, but later we assume this - if (nloc.gt.mnbf) call errquit - & (pname//': nloc > mnbf', 66, UNKNOWN_ERR) - if (nloc.gt.nbf) call errquit - & (pname//': nloc > nbf', 66, UNKNOWN_ERR) - - -c DEBUG some things related to the basis function centers +c DEBUG some things related to the minimal basis centers if (debug) then do a = 1, natoms - if (.not. bas_ce2bfr(basis, a, bflo, bfhi)) - $ call errquit('localized: basis ', 0, BASIS_ERR) - if (master) then - write(luout,'(1x,a,3i4)') 'ao basis, atom = ', a, bflo, bfhi - end if if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) - & call errquit('localized: basis ', 0, BASIS_ERR) + & call errquit(pname//': minbas', 0, BASIS_ERR) if (master) then - write(luout,'(1x,a,3i4)') 'minbas , atom = ', a, bflo, bfhi + write(luout,'(1x,a,3i4)') 'minbas, atom = ', a, bflo, bfhi end if end do end if ! debug -c transform the set of occupied MOs to be localized to the IAO basis -c via C(iao, trans) S C(MO) - - if (.not. ga_create(MT_DBL, mnbf, nloc , 'loc:cibocc', - & mnbf, 0, g_cibocc)) call errquit(pname//': cibocc',0, GA_ERR) - - if (.not. ga_create(MT_DBL, nbf, nloc , 'loc:tmp1', - & nbf, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) - - call ga_dgemm('n', 'n', nbf, nloc, nbf, - & 1.0d0, g_smat, g_c, 0.0d0, g_tmp1) - - call ga_dgemm('t', 'n', mnbf, nloc, nbf, - & 1.0d0, g_iao, g_tmp1, 0.0d0, g_cibocc) - - if (.not. ga_destroy(g_tmp1)) - & call errquit(pname//': ga_destroy failed g_tmp1',0, GA_ERR) - -c if (debug) call ga_print(g_cibocc) - - if (debug .and. master) write(luout,*) 'cibocc done' - -c save a copy of the starting MOs for later - - if (.not.ga_duplicate(g_cibocc, g_tmp2, 'g_tmp2')) call - & errquit(pname//': ga_dup cibocc', 0, GA_ERR) - - call ga_copy(g_cibocc, g_tmp2) - c the MOs in the orthonormal IAO basis must form an orthogonal c matrix. test that before localization, and again later - if (.not. ga_create(MT_DBL, nloc, nloc , 'loc:tmp1', - & mnbf, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) + if (.not. ga_create(MT_DBL, nloc, nloc , 'ibo:tmp', + & 0, 0, g_tmp)) call errquit(pname//': tmp',0, GA_ERR) - call ga_dgemm('t', 'n', nloc, nloc, mnbf, - & 1.0d0, g_cibocc, g_cibocc, 0.0d0, g_tmp1) + call ga_dgemm('t', 'n', nloc, nloc, mnbf, + & 1.0d0, g_c, g_c, 0.0d0, g_tmp) - call ga_add_todiag(g_tmp1, -1.0d0) + call ga_add_todiag(g_tmp, -1.0d0) - call ga_maxelt(g_tmp1, tmp) + call ga_maxelt(g_tmp, tmp) if (master) then - write(luout,'(/1x,a,1x,f18.12)') - & 'IBO loc: largest element of C(MO, trans) C(MO) -1: ', + write(luout,'(/1x,a,1x,f18.8)') + & 'IBO loc: largest element in C(MO,T) C(MO) -1: ', & tmp write(luout,'(1x,a)') & 'should be zero, for CMOs in the IAO basis' end if - if (.not. ga_destroy(g_tmp1)) - & call errquit(pname//': ga_destroy failed g_tmp1',0, GA_ERR) + if (.not. ga_destroy(g_tmp)) + & call errquit(pname//': ga_destroy failed g_tmp',0, GA_ERR) - -c tol = 1d-8 gamma_tol = 1d-10 c if (master) then - write(luout,2) - 2 format(/10x,' iter Max. delocal Mean delocal Converge'/ - $ 10x,' ---- ------------ ------------ ---------') - call util_flush(luout) + write(luout,9001) + 9001 format(/10x,' iter Max. delocal Mean delocal Converge'/ + & 10x,' ---- ------------ ------------ --------') + call util_flush(luout) end if c dprev = 0.0d0 @@ -832,25 +791,24 @@ c c c Analyze convergence by forming P c - d = 0.0d0 - dmax = 0.0d0 - do ss = 1+ga_nodeid(), nloc, ga_nnodes() - s = iloc(ss) - call ga_get(g_cibocc, 1, mnbf, s, s, c(1,1), 1) - qs = 0.0d0 - do a = 1, natoms - if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) - & call errquit(pname//': basis ', 1, BASIS_ERR) - qas = 0.0d0 - do u = bflo, bfhi - qas = qas + c(u,1)*c(u,1) - end do - qs = qs + qas**2 + d = 0.0d0 + dmax = 0.0d0 + do s = 1+ga_nodeid(), nloc, ga_nnodes() + call ga_get(g_c, 1, mnbf, s, s, c(1,1), 1) + qs = 0.0d0 + do a = 1, natoms + if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) + & call errquit(pname//': basis ', 1, BASIS_ERR) + qas = 0.0d0 + do u = bflo, bfhi + qas = qas + c(u,1)*c(u,1) end do -* write(luout,*) ' ds ', s, 1.0d0/qs - dmax = max(dmax, 1.0d0/qs) - d = d + 1.0d0/qs - end do + qs = qs + qas**2 + end do +* write(luout,*) ' ds ', s, 1.0d0/qs + dmax = max(dmax, 1.0d0/qs) + d = d + 1.0d0/qs + end do c #ifdef NWCHEM_USE_GOP_ABSMAX call ga_dgop(1, gamma_max, 1, 'absmax') @@ -861,84 +819,77 @@ c call ga_dgop(1, gamma_max, 1, 'max') call ga_dgop(1, dmax, 1, 'max') #endif - call ga_dgop(2, d , 1, '+') + call ga_dgop(2, d , 1, '+') c - if (ga_nodeid() .eq. 0) then - write(luout,1) iter, dmax, d/dble(nloc), gamma_max - 1 format(10x, i5, 2f15.10, 1p,d12.2,d12.2) - call util_flush(luout) - end if - call ga_sync + if (ga_nodeid() .eq. 0) then + write(luout,'(10x, i5, 2f15.10, 1p,d12.2,d12.2)') + & iter, dmax, d/dble(nloc), gamma_max + call util_flush(luout) + end if + call ga_sync + + if (iter.gt.1 .and. gamma_max.lt.tol) goto 1000 + gamma_max = 0.0d0 + +c Loop over pairs with as much parallelism as possible + + neven = nloc + mod(nloc,2) + do set = 1, neven-1 + do pair = 1+ga_nodeid(), neven/2, ga_nnodes() + call localize_pairs(neven, set, pair, s, t) + if (t .le. nloc) then + + call ga_get(g_c, 1, mnbf, s, s, c(1,1), 1) + call ga_get(g_c, 1, mnbf, t, t, c(1,2), 1) c -** if (abs(d-dprev)/dble(nloc) .lt. tol) goto 1000 -** dprev = d - if (iter.gt.1 .and. gamma_max.lt.tol) goto 1000 - gamma_max = 0.0d0 +c Form rotation information c -c Loop over pairs with as much parallelism as possible + ast = 0.0d0 + bst = 0.0d0 + do a = 1, natoms + if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) + & call errquit(pname//': basis ', 2, + & BASIS_ERR) c - neven = nloc + mod(nloc,2) - do set = 1, neven-1 - do pair = 1+ga_nodeid(), neven/2, ga_nnodes() - call localize_pairs(neven, set, pair, ss, tt) - if (tt .le. nloc) then - s = iloc(ss) - t = iloc(tt) -* write(luout,*) nloc, neven, set, pair, ss, tt, s, t - call ga_get(g_cibocc, 1, mnbf, s, s, c(1,1), 1) - call ga_get(g_cibocc, 1, mnbf, t, t, c(1,2), 1) -c -c Form rotation information -c - ast = 0.0d0 - bst = 0.0d0 - do a = 1, natoms - if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) - & call errquit(pname//': basis ', 2, - & BASIS_ERR) -c - qast = 0.0d0 - qas = 0.0d0 - qat = 0.0d0 - do u = bflo, bfhi - qast = qast + c(u,2)*c(u,1) + c(u,1)*c(u,2) - qas = qas + c(u,1)*c(u,1) - qat = qat + c(u,2)*c(u,2) - end do - qast = qast * 0.5d0 -c -c ast = ast + qast**2 - 0.25d0*(qas - qat)**2 -c bst = bst + qast*(qas - qat) - tmp = -1.0d0*(qas**4 + qat**4) - tmp = tmp + 6.0d0*((qas**2 + qat**2)* qast**2) - tmp = tmp + qas**3 * qat + qas * qat**3 - ast = ast + 0.25d0*tmp - tmp = 4.0d0 * qast*(qas**3 - qat**3) - bst = bst + 0.25d0*tmp - end do -c - gamma = 0.25d0*acos(-ast/sqrt(ast**2+bst**2)) - gamma = sign(gamma,bst) - gamma_max = max(gamma_max, abs(gamma)) -* if (iter .eq. 1 .and. abs(gamma).lt.0.01d0) then -* gamma = (util_random(0)-0.5d0)*3.14d0 -* endif -c - if (abs(gamma) .gt. gamma_tol) then - nrot = nrot + 1 - cosg = cos(gamma) - sing = sin(gamma) -c -c Do the rotation of C -c - call drot(mnbf, c(1,1), 1, c(1,2), 1, cosg, sing) - call ga_put(g_cibocc, 1, mnbf, s, s, c(1,1), 1) - call ga_put(g_cibocc, 1, mnbf, t, t, c(1,2), 1) - end if - end if - end do - call ga_sync - end do + qast = 0.0d0 + qas = 0.0d0 + qat = 0.0d0 + do u = bflo, bfhi + qast = qast + c(u,2)*c(u,1) + c(u,1)*c(u,2) + qas = qas + c(u,1)*c(u,1) + qat = qat + c(u,2)*c(u,2) + end do + qast = qast * 0.5d0 +c +c ast = ast + qast**2 - 0.25d0*(qas - qat)**2 +c bst = bst + qast*(qas - qat) + tmp = -1.0d0*(qas**4 + qat**4) + tmp = tmp + 6.0d0*((qas**2 + qat**2)* qast**2) + tmp = tmp + qas**3 * qat + qas * qat**3 + ast = ast + 0.25d0*tmp + tmp = 4.0d0 * qast*(qas**3 - qat**3) + bst = bst + 0.25d0*tmp + end do +c + gamma = 0.25d0*acos(-ast/sqrt(ast**2+bst**2)) + gamma = sign(gamma,bst) + gamma_max = max(gamma_max, abs(gamma)) + + if (abs(gamma) .gt. gamma_tol) then + nrot = nrot + 1 + cosg = cos(gamma) + sing = sin(gamma) +c +c Do the rotation of C +c + call drot(mnbf, c(1,1), 1, c(1,2), 1, cosg, sing) + call ga_put(g_c, 1, mnbf, s, s, c(1,1), 1) + call ga_put(g_c, 1, mnbf, t, t, c(1,2), 1) + end if + end if + end do + call ga_sync + end do end do c 1000 continue @@ -947,178 +898,29 @@ c c the MOs in the orthonormal IAO basis must form an orthogonal c matrix. test that after localization - if (.not. ga_create(MT_DBL, nloc, nloc , 'loc:tmp1', - & mnbf, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) + if (.not. ga_create(MT_DBL, nloc, nloc , 'ibo:tmp', + & mnbf, 0, g_tmp)) call errquit(pname//': tmp',0, GA_ERR) - call ga_dgemm('t', 'n', nloc, nloc, mnbf, - & 1.0d0, g_cibocc, g_cibocc, 0.0d0, g_tmp1) + call ga_dgemm('t', 'n', nloc, nloc, mnbf, + & 1.0d0, g_c, g_c, 0.0d0, g_tmp) - call ga_add_todiag(g_tmp1, -1.0d0) + call ga_add_todiag(g_tmp, -1.0d0) - call ga_maxelt(g_tmp1, tmp) + call ga_maxelt(g_tmp, tmp) if (master) then - write(luout,'(/1x,a,1x,f18.12)') - & 'IBO loc: largest element of C(MO, trans) C(MO) -1: ', + write(luout,'(/1x,a,1x,f18.8)') + & 'IBO loc: largest element in C(MO,T) C(MO) -1: ', & tmp write(luout,'(1x,a)') - & 'should be zero, for LMOs in the IAO basis' + & 'should be zero, for IBOs in the IAO basis' end if -c Following steps: - -c 1. while we have array g_tmp1 allocated, also calculate the -c localization transform. The starting MOs were saved in g_tmp2 in -c the IAO basis + if (.not. ga_destroy(g_tmp)) + & call errquit(pname//': ga_destroy failed g_tmp',0, GA_ERR) -c 2. calculate LMO energies and order the LMO set by energy + return + end -c 3. save the localization transform to a file, instead of -c generating it again in the calling routine and saving it from -c there - -c 4. We then transform the CMOs to re-ordered LMOs in the -c AO basis and store the transformed MOs back in g_c - - -c CMOs(iao, trans) * LMOs(iao) = localization transform -> g_tmp1 - call ga_dgemm('t', 'n', nloc, nloc, mnbf, - & 1.0d0, g_tmp2, g_cibocc, 0.0d0, g_tmp1) - - if (.not. ga_destroy(g_tmp2)) - & call errquit(pname//': ga_destroy failed g_tmp2',0, GA_ERR) - -c transform MO energies. We assume that the first nloc MOs were -c localized, and that the Fock matrix was diagonal in the basis of -c input MOs. We also assume (this is tested upon entry) that -c nloc is smaller than the column sizes of arays eval, c(:,1:2) -c so we don't segfault when using c(:,:) as temp memory - - do jmo = 1,nloc - call ga_get(g_tmp1, 1, nloc, jmo, jmo, c(1,1), 1) - tmp = 0.0d0 - do imo = 1,nloc - tmp = tmp + eval(imo) * c(imo,1)**2 - end do ! jmo - c(jmo,2) = tmp ! temp storage for transformed MO energies - end do ! imo - do imo = 1,nloc - eval(imo) = c(imo,2) - end do - -c sort LMOs by energy. We'll have to do this on the master node, or -c otherwise avoid having multiple nodes working on the same data - - if (master) then - do imo = 1,nloc -1 - minval = eval(imo) -c find lowest eval(kmo) below eval(imo) - kmo = 0 - do jmo = imo+1, nloc - if (eval(jmo).lt.minval) then - kmo = jmo - minval = eval(jmo) - end if - end do -c if kmo > 0 we swap imo and jmo - if (kmo.ne.0) then - if (debug.and.master) write(luout,*) 'swap ',imo,kmo -c swap eigenvalues - swap = eval(imo) - eval(imo) = eval(kmo) - eval(kmo) = swap -c swap columns of LMO transformation - call ga_get(g_tmp1, 1, nloc, imo, imo, c(1,1), 1) - call ga_get(g_tmp1, 1, nloc, kmo, kmo, c(1,2), 1) - call ga_put(g_tmp1, 1, nloc, imo, imo, c(1,2), 1) - call ga_put(g_tmp1, 1, nloc, kmo, kmo, c(1,1), 1) -c swap LMOs in IAO basis, for the analysis below - call ga_get(g_cibocc, 1, mnbf, imo, imo, c(1,1), 1) - call ga_get(g_cibocc, 1, mnbf, kmo, kmo, c(1,2), 1) - call ga_put(g_cibocc, 1, mnbf, imo, imo, c(1,2), 1) - call ga_put(g_cibocc, 1, mnbf, kmo, kmo, c(1,1), 1) - end if - end do - end if ! master - - call ga_sync - -c write transformation to scratch file - - call util_file_name('lmotrans',.true.,.true.,lmotrans) - if(.not.file_write_ga(lmotrans,g_tmp1)) call errquit - & (pname//': could not write lmotrans',0, DISK_ERR) - - if (debug.and.master) write (luout,*) - & pname//': lmos reordered and lmotrans written' - -c transform CMOs to re-ordered LMOs - - if (.not. ga_create(MT_DBL, nbf, nloc , 'loc:tmp2', - & nbf, 0, g_tmp2)) call errquit(pname//': tmp2',0, GA_ERR) - - call ga_dgemm('n', 'n', nbf, nloc, nloc, - & 1.0d0, g_c, g_tmp1, 0.0d0, g_tmp2) - - call ga_copy(g_tmp2, g_c) - - if (debug.and.master) write (luout,*) pname//': g_c now has LMOs' - - if (.not. ga_destroy(g_tmp1)) - & call errquit(pname//': ga_destroy failed g_tmp1',0, GA_ERR) - - if (.not. ga_destroy(g_tmp2)) - & call errquit(pname//': ga_destroy failed g_tmp2',0, GA_ERR) - -c Analyze localization of each MO: -c per LMO, a list of atomic populations is printed -c in decreasing magnitude, with the polulations in parentheses - - if (master) then - write(luout,'(/1x,a/)') 'IAO-IBO localized orbitals' - do ss = 1, nloc - s = iloc(ss) - call ga_get(g_cibocc, 1, mnbf, s, s, c(1,1), 1) - nlist = 0 - do a = 1, natoms - if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) - & call errquit('localized: basis ', 3, - & BASIS_ERR) - qas = 0.0d0 - do u = bflo, bfhi - qas = qas + c(u,1)*c(u,1) - end do - if (abs(qas) .gt. 0.01d0) then - nlist = nlist + 1 - list(nlist) = a - pop(nlist) = qas - end if - end do - do u = 1, nlist - do t = 1, u-1 - if (abs(pop(t)).lt.abs(pop(u))) then - tmp = pop(u) - pop(u) = pop(t) - pop(t) = tmp - tt = list(u) - list(u) = list(t) - list(t) = tt - end if - end do - end do - write(luout,77) s, eval(s), - & occ(s),(list(a), pop(a), a=1,nlist) - 77 format(i5, 1x, f14.6,1x, f5.3, 1x,100(2x,i4,'(',f5.2,')')) - end do - write(luout,*) - call util_flush(luout) - end if - - if (.not. ga_destroy(g_cibocc)) - & call errquit(pname//': ga_destroy failed g_cibocc',0, GA_ERR) - - if (debug.and.master) write(luout,*) 'leaving '//pname - - end c c ================================================================= c diff --git a/src/property/GNUmakefile b/src/property/GNUmakefile index e85f2e9bfb..9a0611c846 100644 --- a/src/property/GNUmakefile +++ b/src/property/GNUmakefile @@ -63,6 +63,7 @@ CalcPerturbedTDPmat1.o \ CalcPerturbedTDPmat1_opt.o \ localization_driver.o \ + ibo_localization.o \ hnd_vec_write.o \ giao_b1_movecs_tools.o \ aor_r1_beta_anl_tools.o\ diff --git a/src/property/ibo_localization.F b/src/property/ibo_localization.F new file mode 100644 index 0000000000..db5b52c44b --- /dev/null +++ b/src/property/ibo_localization.F @@ -0,0 +1,967 @@ + + subroutine ibo_localization(rtdb, geom, ltyp, basis, g_movecs, + & nocc, nvir, nmo, nbf, natoms, eval, occ, c, pop, list) + +c ================================================================= +c IAO construction and generation of occupied or virtual IBOs. +c Written by J. Autschbach, jochena@buffalo.edu +c +c The formalism is described in G. Knizia, J. Chem. Theory +c Comput. 2013, 9, 4834 (with additional useful information provided +c in Senjean et al., J. Chem. Theory Comput. 2021, 17, 1337 and the +c SI of that article). What is implemented here corresponds to +c Knizia's original formulation, and the equations for p = 4 in +c Appendix D as given in the 'bugfixed' version of the article on +c the author's web page [https://sites.psu.edu/knizia/software/ +c (accessed 02/2023)]. + +c g_movecs and arrays occ and eval are for the MOs of a given +c spin (alpha or beta) + +c The localization routine called from here is in +c src/ddscf/localize.F + +c ================================================================= + + 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" +#include "apiP.fh" +#include "prop.fh" +#include "bgj.fh" + + +c subroutine arguments: + + integer rtdb, geom, basis + character*(3) ltyp + integer g_movecs + integer nocc, nvir, nmo, nbf, natoms + double precision eval(nbf), occ(nbf), c(nbf,2) + double precision pop(natoms) + integer list(natoms) + +c local GA handles: + + integer g_s2, g_s12, g_p12, g_p21 + integer g_ctilde, g_iao, g_mo + integer g_temp, g_tmp1, g_tmp2, g_cib, g_u, g_vt + integer g_smat + +c local variables: + + logical master, debug + integer dbg, info + + integer ncore, imo, jmo, kmo + + integer n1, n2, nsing, noff + + integer s, nlist, a, i, t, tt, u, bflo, bfhi + + double precision rtemp, qas + + double precision small + parameter (small=1d-8) + + integer + & mnbf, iao_mxprim, iao_mxang, iao_mxcont, + & iao_mxnbf_cn, iao_mxnbf_ce, iao_nshells, + & minbas + + double precision minval, swap + integer l_val, k_val + + logical int_normalize + external int_normalize + + integer ga_create_atom_blocked + external ga_create_atom_blocked + + + character*(256) lmotrans + logical file_write_ga + external file_write_ga + + character*(16) pname + +c ================================================================= + + pname = 'ibo_localization' + + dbg = 0 + master = ga_nodeid().eq.0 + debug = (dbg>0) .and. master ! .true. during development + + if (ltyp.ne.'occ' .and. ltyp.ne.'vir') call errquit + & (pname//': loc. type unknown', 0, BASIS_ERR) + + if(debug) then + if (ltyp.eq.'occ') write (luout,*) + & 'entering occupied IBO localization' + if (ltyp.eq.'vir') write (luout,*) + & 'entering virtual IBO localization' + end if + + if (.not. geom_num_core(rtdb, geom, 'ddscf', ncore)) ncore = 0 + + if (debug) write (luout,*) 'ncore = ',ncore + + + +c ------------------------------ +c construct IAOs in the AO basis +c ------------------------------ + +c Note: +c Basis 1 is the AO basis used in the SCF calculation +c Basis 2 is the minimal auxiliary basis + +c AO Overlap Matrix S1 -> g_smat: + + 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 auxiliary basis overlap S2 -> g_s2 basis needs to defined in +c the input as "iao basis". we will now create the basis here and +c then calculate the overlap and the mixed ao-iao basis overlap +c S12 + + if (.not. bas_create(minbas, 'iao basis')) + & call errquit(pname//': cannot create iao bas', 86, BASIS_ERR) + + if (.not. bas_rtdb_load(rtdb, geom, minbas, 'iao basis')) + & call errquit(pname//': iao basis not present', 86, BASIS_ERR) + + if (.not. int_normalize(rtdb,minbas)) + & call errquit(pname//': int_normalize iao', 86, INT_ERR) + + if (.not. bas_numbf(minbas, mnbf)) + & call errquit(pname//': bas err iao mnbf', 86, BASIS_ERR) + + if (.not. bas_nprim_cn_max(minbas,iao_mxprim)) + & call errquit(pname//': bas err iao mxprim', 86, BASIS_ERR) + + if (.not. bas_high_angular(minbas,iao_mxang)) + & call errquit(pname//': bas err iao mxang', 86, BASIS_ERR) + + if (.not. bas_ncontr_cn_max(minbas,iao_mxcont)) + & call errquit(pname//': bas err iao mxcont', 86, BASIS_ERR) + + if (.not. bas_nbf_cn_max(minbas,iao_mxnbf_cn)) + & call errquit(pname//': bas err iao mcnbf_cn', 86, BASIS_ERR) + + if (.not. bas_nbf_ce_max(minbas,iao_mxnbf_ce)) + & call errquit(pname//': bas err iao mxnbf_ce', 86, BASIS_ERR) + + if (.not. bas_numcont(minbas,iao_nshells)) + & call errquit(pname//': bas err iao nshells', 86, BASIS_ERR) + +c not sure if it makes much sense in case the minimal basis for +c the IAOs is larger than our regular basis, so let's check that +c here and exit with an error if nbf < mnbf + + if (nbf.lt.mnbf) call errquit + & (pname//': nbf < mnbf. cannot handle!', 66, UNKNOWN_ERR) + +c create overlap for minbas, and the mixed basis-minbas +c overlap S12. +c Then we calculate the projectors P12 and P21 +c via solving S1 P12 = S12 for P12 etc. + + g_s2 = ga_create_atom_blocked(geom, minbas, 'loc:s2') + call ga_zero(g_s2) + call int_1e_ga(minbas, minbas, g_s2, 'overlap', .false.) + + if (.not. ga_create(MT_DBL, nbf, mnbf, 'loc:s12', + & nbf, 0, g_s12)) call errquit(pname//': S12',61, GA_ERR) + + call ga_zero(g_s12) + call int_1e_ga(basis, minbas, g_s12, 'overlap', .false.) + + if (.not.ga_duplicate(g_s12, g_p12, 'g_p12')) call + & errquit(pname//': ga_dup P12',62,GA_ERR) + + call ga_copy(g_s12, g_p12) + +#ifdef SCALAPACK + info = ga_solve(g_smat, g_p12) + if (debug) write(luout,*) + & 'info from ga_solve P12:', info +#else + if (debug) write(luout,*) 'call 1 ma_solve' + call ma_solve(g_smat, g_p12) +#endif + + if (.not. ga_create(MT_DBL, mnbf, nbf, 'loc:p21', + & mnbf, 0, g_p21)) call errquit(pname//': P21',63, GA_ERR) + + call ga_transpose(g_s12, g_p21) + if (debug) write(luout,*) 'P12 assembled' + +c S12 no longer needed + if (.not. ga_destroy(g_s12)) + & call errquit(pname//': ga_destroy failed g_s12',61, GA_ERR) + +c : scalapack throws a 'ga_to_SL: rows/cols error' so we're using +c ma_solve for the time being +c#ifdef SCALAPACK +c info = ga_solve(g_s2, g_p21) +c if (debug) write(luout,*) +c & 'info from ga_solve P21:', info +c#else + if (debug) write(luout,*) 'call 2 ma_solve' + call ma_solve(g_s2, g_p21) +c#endif + +c g_s2 no longer needed + if (.not. ga_destroy(g_s2)) call errquit( + & pname//': error destroying g_s2',0, GA_ERR) + + if (dbg>1) then + call ga_print(g_p12) + call ga_print(g_p21) + end if + +c store P12 * P21 in array g_temp; we'll need it soon + + if (.not.ga_duplicate(g_smat, g_temp, 'g_temp')) call + & errquit(pname//': ga_dup temp', 71, GA_ERR) + + call ga_dgemm('n', 'n', nbf, nbf, mnbf, + & 1.0d0, g_p12, g_p21, 0.0d0, g_temp) + + if (debug) write (luout,*) 'p12 p21 -> temp' + +c P21 is no longer needed + + if (.not. ga_destroy(g_p21)) + & call errquit(pname//': ga_destroy failed g_p21',61, GA_ERR) + +c construct IAOS from occ. MOs: + +c copy the relevant CMOs to g_mo + + if (.not. ga_create(MT_DBL, nbf, nocc, 'loc:g_mo', + & nbf, 0, g_mo)) call errquit(pname//': g_mo',0, GA_ERR) + + call ga_copy_patch('n', + & g_movecs, 1, nbf, 1, nocc, + & g_mo, 1, nbf, 1, nocc) + + if (debug) write (luout,*) 'movecs(occ) -> mo' + +c create C-tilde from Appendix C of Knizia's IBO paper. g_temp +c holds P12 * P21; we won't need it after the next matrix +c multiplication + + if (.not. ga_create(MT_DBL, nbf, nocc, 'loc:ctilde', + & nbf, 0, g_ctilde)) call errquit('loc_driver: sc',0, GA_ERR) + + call ga_dgemm('n', 'n', nbf, nocc, nbf, + & 1.0d0, g_temp, g_mo, 0.0d0, g_ctilde) + + if (.not. ga_destroy(g_temp)) + & call errquit(pname//': ga_destroy failed g_temp',71, GA_ERR) + + if (debug) write (luout,*) 'pre sym_orth 1' + call sym_orth(g_ctilde, g_smat, nbf, nocc) + + if (dbg>1) call ga_print(g_ctilde) + +c now assemble the matrix products to form matrix A from Appendix +c C of the IBO paper. The result will be accumulated in array +c g_iao which are the IAOs in basis 1. The shape of the array is +c the same as P12 + + if (.not.ga_duplicate(g_smat, g_tmp1, 'g_tmp1')) call + & errquit(pname//': ga_dup tmp1', 72, GA_ERR) + + if (.not.ga_duplicate(g_smat, g_tmp2, 'g_tmp2')) call + & errquit(pname//': ga_dup tmp2', 72, GA_ERR) + + if (.not.ga_duplicate(g_smat, g_temp, 'g_temp')) call + & errquit(pname//': ga_dup temp', 72, GA_ERR) + + if (.not.ga_duplicate(g_p12, g_iao, 'g_iao')) call + & errquit(pname//': ga_dup iao', 72, GA_ERR) + + +c C * C(T) + call ga_dgemm('n', 't', nbf, nbf, nocc, + & 1.0d0, g_mo, g_mo, 0.0d0, g_temp) + +c C * C(T) * S1 -> tmp1 + call ga_dgemm('n', 'n', nbf, nbf, nbf, + & 1.0d0, g_temp, g_smat, 0.0d0, g_tmp1) + +c tilde-C * tilde-C(T) + call ga_dgemm('n', 't', nbf, nbf, nocc, + & 1.0d0, g_ctilde, g_ctilde, 0.0d0, g_temp) + +c tilde-C * tilde-C(T) * S1 -> tmp2 + call ga_dgemm('n', 'n', nbf, nbf, nbf, + & 1.0d0, g_temp, g_smat, 0.0d0, g_tmp2) + +c C * C(T) * S1 * tilde-C * tilde-C(T) * S1 + + call ga_dgemm('n', 'n', nbf, nbf, nbf, + & 1.0d0, g_tmp1, g_tmp2, 0.0d0, g_temp) + +c multiply prev. result by P12 and store in g_iao + + call ga_dgemm('n', 'n', nbf, mnbf, nbf, + & 1.0d0, g_temp, g_p12, 0.0d0, g_iao) + +c form C * C(T) -1 + call ga_add_todiag(g_tmp1, -1.0d0) + +c form tilde-C * tilde-C(T) -1 + call ga_add_todiag(g_tmp2, -1.0d0) + +c multiply the two prev. results with each other, then multiply +c the result with P12 and add to g_iao + + call ga_dgemm('n', 'n', nbf, nbf, nbf, + & 1.0d0, g_tmp1, g_tmp2, 0.0d0, g_temp) + + call ga_dgemm('n', 'n', nbf, mnbf, nbf, + & 1.0d0, g_temp, g_p12, 1.0d0, g_iao) + +c g_ctilde no longer needed + + if (.not. ga_destroy(g_ctilde)) + & call errquit(pname//': ga_destroy failed ctilde',62, GA_ERR) + +c delete temp arrays before calling sym_orth, which allocates its +c own temp memory + + if (.not. ga_destroy(g_tmp1)) + & call errquit(pname//': ga_destroy failed g_tmp1',71, GA_ERR) + + if (.not. ga_destroy(g_tmp2)) + & call errquit(pname//': ga_destroy failed g_tmp2',71, GA_ERR) + + if (.not. ga_destroy(g_temp)) + & call errquit(pname//': ga_destroy failed g_temp',71, GA_ERR) +c last step: orthonormalize g_iao + + if (debug) write (luout,*) 'pre sym_orth 2' + call sym_orth(g_iao, g_smat, nbf, mnbf) + + call ga_sync + + if (dbg>1) call ga_print(g_iao) + +c check orthogonality of IAOs w.r.t. metric S1 +c g_p12 is used for scratch + + if (.not. ga_create(MT_DBL, mnbf, mnbf, 'g_temp', + & mnbf, 0, g_temp)) call errquit(pname//': temp',73, GA_ERR) + + call ga_dgemm('n', 'n', nbf, mnbf, nbf, + & 1.0d0, g_smat, g_iao, 0.0d0, g_p12) + + call ga_dgemm('t', 'n', mnbf, mnbf, nbf, + & 1.0d0, g_iao, g_p12, 0.0d0, g_temp) + + if (dbg>1) call ga_print(g_temp) + + call ga_add_todiag(g_temp, -1.0d0) + + call ga_maxelt(g_temp, rtemp) + if (master) then + write(luout,'(/1x,a,1x,f18.8)') + & 'IBO loc: largest element in C(iao,T) S C(iao) -1: ', + & rtemp + write(luout,'(1x,a/1x,a/)') + & 'Significant deviations from zero may indicate', + & 'elevated numerical noise in the IAO generation' + end if + + if (.not. ga_destroy(g_temp)) + & call errquit(pname//': ga_destroy failed g_temp',73, GA_ERR) + +c g_p12 no longer needed + + if (.not. ga_destroy(g_p12)) + & call errquit(pname//': ga_destroy failed g_p12',61, GA_ERR) + +c --------------------------------------------------------------- +c IAOs are now in array g_iao. Next, generate occupied or virtual +c IBOs, depending on the input settings (ltyp) +c --------------------------------------------------------------- + + if (ltyp.eq.'occ') then + +c ----------------------- +c generate occupied IBOs: +c ----------------------- + +c note: g_mo already holds the occupied MOs + + if (debug) then + write(luout,*) 'movecs(occ) -> mo' + end if + +c few more sanity check, just in case + if (nocc.gt.mnbf) call errquit + & (pname//': nocc > mnbf', 66, UNKNOWN_ERR) + if (nocc.gt.nbf) call errquit + & (pname//': nocc > nbf', 66, UNKNOWN_ERR) + +c transform the set of occupied MOs to be localized to the IAO +c basis via C(iao,T) S C(MO) + + if (.not. ga_create(MT_DBL, mnbf, nocc , 'loc:cib', + & mnbf, 0, g_cib)) + & call errquit(pname//': cib',0, GA_ERR) + + if (.not. ga_create(MT_DBL, nbf, nocc , 'loc:tmp1', + & 0, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) + + call ga_dgemm('n', 'n', nbf, nocc, nbf, + & 1.0d0, g_smat, g_mo, 0.0d0, g_tmp1) + + call ga_dgemm('t', 'n', mnbf, nocc, nbf, + & 1.0d0, g_iao, g_tmp1, 0.0d0, g_cib) + + if (.not. ga_destroy(g_tmp1)) + & call errquit(pname//': ga_destroy failed g_tmp1',0, GA_ERR) + + if (dbg>1) call ga_print(g_cib) + + if (debug) write(luout,*) 'cib done' + +c save a copy of the starting MOs for later + + if (.not.ga_duplicate(g_cib, g_tmp2, 'g_tmp2')) call + & errquit(pname//': ga_dup cib', 0, GA_ERR) + + call ga_copy(g_cib, g_tmp2) + +c perform localization of the MOs in IAO basis: + + call localizeIBO(minbas, c, g_cib, nocc, nbf, mnbf, + & natoms) + + + else if (ltyp.eq.'vir') then + +c ---------------------- +c generate virtual IBOs: +c ---------------------- + +c store the virtual MOs in g_mo. The array is allocated, so we +c destroy it frst + + if (.not. ga_destroy(g_mo)) call errquit( + & pname//': error destroying g_mo',1, GA_ERR) + + if (.not. ga_create(MT_DBL, nbf, nvir, 'loc:g_mo', + & nbf, 0, g_mo)) call errquit(pname//': g_mo',1, GA_ERR) + + call ga_copy_patch('n', + & g_movecs, 1, nbf, nocc+1, nmo, + & g_mo, 1, nbf, 1, nvir) + +c sanity check: for virtual localization, assume no. of +c virtuals >= (mnbf-nocc) + + if (nvir.lt.(mnbf-nocc)) call errquit + & (pname//': nvir < mnbf-nocc. Cannot handle', 66, + & UNKNOWN_ERR) + + if (debug) then + write(luout,*) 'movecs(vir) -> mo' + write(luout,*) 'nocc, nvir, nmo, nbf, mnbf', + & nocc, nvir, nmo, nbf, mnbf + end if + +c transform the set of virtual MOs to the IAO basis +c via C(iao,T) S C(MO), store in g_tmp2 (and keep the array) + + if (.not. ga_create(MT_DBL, mnbf, nvir , 'loc:tmp2', + & 0, 0, g_tmp2)) call errquit(pname//': tmp2',0, GA_ERR) + + if (.not. ga_create(MT_DBL, nbf, nvir , 'loc:tmp1', + & 0, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) + + call ga_dgemm('n', 'n', nbf, nvir, nbf, + & 1.0d0, g_smat, g_mo, 0.0d0, g_tmp1) + + call ga_dgemm('t', 'n', mnbf, nvir, nbf, + & 1.0d0, g_iao, g_tmp1, 0.0d0, g_tmp2) + + if (.not. ga_destroy(g_tmp1)) + & call errquit(pname//': ga_destroy failed g_tmp1',0, GA_ERR) + + if (debug .and. master) write(luout,*) 'cib before svd' + +c now perform SVD on the matrix C(vir) that we just calculated and +c stored in array g_tmp2. The "eigenvector" matrix U gives the +c true virtual MOs in the IAO basis that we want to localize + +c #(sing. values) = min(mnbf,nvir) = mnbf , assuming virtuals. +c The non-zero singular values must number mnbf - nocc, which +c means we will localize a number of virtuals that's smaller than +c nvir and also smaller than min(mnbf,nvir). + + nsing=min(mnbf,nvir) + + if (.not. ma_push_get(MT_DBL, nsing, 'svd-vals', l_val, k_val)) + & call errquit(pname//': svd-vals', nsing, MA_ERR) + + if (.not. ga_create(MT_DBL, mnbf, nsing,'g_u', 0, 0, g_u)) + & call errquit(pname//': svd-U', mnbf*nsing, GA_ERR) + + if (.not. ga_create(MT_DBL, nsing, nvir, 'g_vt',0, 0, g_vt)) + & call errquit(pname//': svd-vT', nsing*nvir, GA_ERR) + + if (debug) write(luout,*) 'svd: m,n,nsing=',mnbf,nvir,nsing + + call wrap_dgesvd(g_tmp2, g_u, g_vt, dbl_mb(k_val)) + +c wrap_dgesvd leaves the content of the array in arg 1 intact, no +c matter what the options in the call to dgesvd are inside the +c routine, so g_tmp2 still contains the initial (redundant) +c virtual MO set expressed in the IAO basis. We'll use it later. + + if (debug) then + write(luout,*) pname//' singular values:' + do imo = 1,mnbf + write(luout,*) imo, dbl_mb(k_val+imo-1) + end do + end if + +c count the non-zero singular values: re-define nsing + + nsing = 0 + do imo = 1,mnbf + if (dbl_mb(k_val+imo-1) .gt.small) nsing = nsing+1 + end do + if (master) write(luout,'(1x,a,i10)') + & 'non-zero singular values: ',nsing + +c we should have as many non-zero singlular values as we have +c additional IAOs beyond the number of occupied MOs. Check that, +c and print a warning if the values differ. + + if (nsing .ne. (mnbf-nocc)) then + if (master) + & write(luout,'(/1x,a/1x,a/1x,a/1x,a,i6/1x,a,i6/1x,a)') + & 'WARNING: number of non-zero singular values differs from', + & ' the number of IAOs minus the number of occupied MOs.', + & ' Will continue suspiciously!', + & ' nsing : ', nsing, + & ' mnbf - nocc : ', mnbf - nocc, + & ' The numbers in the two previous lines should be the same.' + end if + +c array g_cib will hold the left-singular vectors corresponding to +c non-zero singlular values. Since dgesvd returns them in +c descending order, and they are by definition positive, we simply +c take the first nsing vectors + + if (.not. ga_create(MT_DBL, mnbf, nsing,'loc:cib', 0, 0, g_cib)) + & call errquit(pname//': cib', mnbf*nsing, GA_ERR) + + call ga_copy_patch('n', + & g_u, 1, mnbf, 1, nsing, + & g_cib, 1, mnbf, 1, nsing ) + +c dellocate memory used for SVD + + if (.not. ga_destroy(g_vt)) + & call errquit(pname//': ga_destroy failed g_vt',0, GA_ERR) + + if (.not. ga_destroy(g_u)) + & call errquit(pname//': ga_destroy failed g_u',0, GA_ERR) + + if (.not. ma_pop_stack(l_val)) call errquit( + & pname//': error MA pop svd-vals',0, MA_ERR) + + if (dbg>1) call ga_print(g_cib) + + if (debug) write(luout,*) 'cib done' + +c perform localization of the MOs in IAO basis: + + call localizeIBO(minbas, c, g_cib, nsing, nbf, mnbf, + & natoms) + + end if ! ltyp occ or vir + +c localization interations done + +c calculate the localization transform. The starting MOs were +c saved in g_tmp2 in the IAO basis +c CMOs(iao,T) * LMOs(iao) = localization transform -> g_tmp1 + + if (ltyp.eq.'occ') then + + if (.not. ga_create(MT_DBL, nocc, nocc, 'loc:tmp1', + & 0, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) + + n1 = nocc ! dimensioning used later + n2 = nocc + noff = 0 + + call ga_dgemm('t', 'n', nocc, nocc, mnbf, + & 1.0d0, g_tmp2, g_cib, 0.0d0, g_tmp1) + + if (debug) write(luout,*) 'ibo occ: transformation in g_tmp1' + + elseif (ltyp.eq.'vir') then + + n1 = nvir ! dimensioning used later + n2 = nsing + noff = nocc ! offset for virtual orbitals + + if (.not. ga_create(MT_DBL, nvir, nsing, 'loc:tmp1', + & 0, 0, g_tmp1)) call errquit(pname//': tmp1',0, GA_ERR) + + call ga_dgemm('t', 'n', nvir, nsing, mnbf, + & 1.0d0, g_tmp2, g_cib, 0.0d0, g_tmp1) + + end if ! ltyp + + if(noff+n1 > nbf) call errquit (pname// + & ': noff + n1 > nbf', noff, UNKNOWN_ERR) + + if (.not. ga_destroy(g_tmp2)) + & call errquit(pname//': ga_destroy failed g_tmp2',0, GA_ERR) + +c transform MO energies. We assume that the Fock matrix was diagonal +c in the basis of input MOs. We also assume that the loop dimensions +c do not exceed the capacity of array c(:,:) + + do jmo = 1,n2 + call ga_get(g_tmp1, 1, n1, jmo, jmo, c(1,1), 1) + rtemp = 0.0d0 + do imo = 1,n1 + rtemp = rtemp + eval(imo+noff) * c(imo,1)**2 + end do ! jmo + c(jmo,2) = rtemp ! transformed MO energies + end do ! imo + do imo = 1,n2 + eval(imo+noff) = c(imo,2) + end do + + if (debug) write(luout,*) 'transformed MO energies' + +c sort IBOs by energy. We'll have to do this on the master +c node, or otherwise avoid having multiple nodes working on the +c same data + + if (master) then + do imo = 1,n2 -1 + minval = eval(imo+noff) +c find lowest eval(kmo) below eval(imo) (with offsets) + kmo = 0 + do jmo = imo+1, n2 + if (eval(jmo+noff).lt.minval) then + kmo = jmo + minval = eval(jmo+noff) + end if + end do +c if kmo > 0 we swap imo and jmo + if (kmo.ne.0) then + if (debug) write(luout,*) 'swap ',imo,kmo +c swap eigenvalues + swap = eval(imo+noff) + eval(imo+noff) = eval(kmo+noff) + eval(kmo+noff) = swap +c swap columns of LMO transformation + call ga_get(g_tmp1, 1, n1, imo, imo, c(1,1), 1) + call ga_get(g_tmp1, 1, n1, kmo, kmo, c(1,2), 1) + call ga_put(g_tmp1, 1, n1, imo, imo, c(1,2), 1) + call ga_put(g_tmp1, 1, n1, kmo, kmo, c(1,1), 1) +c swap LMOs in IAO basis, for the analysis below + call ga_get(g_cib, 1, mnbf, imo, imo, c(1,1), 1) + call ga_get(g_cib, 1, mnbf, kmo, kmo, c(1,2), 1) + call ga_put(g_cib, 1, mnbf, imo, imo, c(1,2), 1) + call ga_put(g_cib, 1, mnbf, kmo, kmo, c(1,1), 1) + end if + end do + end if ! master + + call ga_sync + +c write transformation to scratch file + + call util_file_name('lmotrans',.true.,.true.,lmotrans) + if(.not.file_write_ga(lmotrans,g_tmp1)) call errquit + & (pname//': could not write lmotrans',0, DISK_ERR) + + if (debug) write (luout,*) + & pname//': lmos reordered and lmotrans written' + +c transform CMOs to re-ordered LMOs + + if (.not. ga_create(MT_DBL, nbf, n2 , 'loc:tmp2', + & nbf, 0, g_tmp2)) call errquit(pname//': tmp2',0, GA_ERR) + + call ga_dgemm('n', 'n', nbf, n2, n1, + & 1.0d0, g_mo, g_tmp1, 0.0d0, g_tmp2) + + + if (debug) write (luout,*) pname//': g_tmp2 now has LMOs' + + if (.not. ga_destroy(g_tmp1)) + & call errquit(pname//': ga_destroy failed g_tmp1',0, GA_ERR) + +c copy the IBOs into the relevant part of movecs + + if (ltyp.eq.'occ') then + call ga_copy_patch('n', + & g_tmp2, 1, nbf, 1, nocc, + & g_movecs, 1, nbf, 1, nocc) + + else if (ltyp.eq.'vir') then + call ga_copy_patch('n', + & g_tmp2, 1, nbf, 1, nsing, + & g_movecs, 1, nbf, nocc+1, nocc+nsing) + end if + + if (debug) write (luout,*) pname//': IBOs copied to g_movecs' + + if (.not. ga_destroy(g_tmp2)) + & call errquit(pname//': ga_destroy failed g_tmp2',0, GA_ERR) + +c Analyze localization of each MO: per LMO, a list of atomic +c populations is printed in decreasing magnitude, with the +c polulations in parentheses. This code is equivalent to the on +c found in the Pipek-Mezey localization routine + + if (master) then + write(luout,'(/1x,a/)') 'IAO-IBO localized orbitals' + do s = 1, n2 + call ga_get(g_cib, 1, mnbf, s, s, c(1,1), 1) + nlist = 0 + do a = 1, natoms + if (.not. bas_ce2bfr(minbas, a, bflo, bfhi)) + & call errquit('localized: basis ', 3, + & BASIS_ERR) + qas = 0.0d0 + do u = bflo, bfhi + qas = qas + c(u,1)*c(u,1) + end do + if (abs(qas) .gt. 0.01d0) then + nlist = nlist + 1 + list(nlist) = a + pop(nlist) = qas + end if + end do + do u = 1, nlist + do t = 1, u-1 + if (abs(pop(t)).lt.abs(pop(u))) then + rtemp = pop(u) + pop(u) = pop(t) + pop(t) = rtemp + tt = list(u) + list(u) = list(t) + list(t) = tt + end if + end do + end do + write(luout,9002) s, eval(s+noff), + & occ(s+noff),(list(a), pop(a), a=1,nlist) + 9002 format(i5, 1x, f14.6,1x, f5.3, 1x,100(2x,i4,'(',f5.2,')')) + end do + write(luout,*) + call util_flush(luout) + end if + + if (.not. ga_destroy(g_cib)) + & call errquit(pname//': ga_destroy failed g_cib',0, GA_ERR) + +c deallocate remaining arrays that are no longer needed + + if (.not. ga_destroy(g_iao)) call errquit( + & pname//': error destroying g_iao',0, GA_ERR) + + if (.not. ga_destroy(g_mo)) call errquit( + & pname//': error destroying g_mo',0, GA_ERR) + +c smat not needed anymore + if (.not. ga_destroy(g_smat)) call errquit( + & pname//': error destroying g_smat',0, GA_ERR) + +c destroy minimal basis (iao basis) + if (.not.bas_destroy(minbas)) + & call errquit(pname//'iao bas_destroy failed',0,BASIS_ERR) + + if (ltyp.eq.'occ') then + n1 = 1 + n2 = nocc + else if (ltyp.eq.'vir') then + n1 = nocc + 1 + n2 = nocc + nsing + end if + + if (master) write(luout, + & '(/1x,a,a,a/1x,a/1x,i10,2x,a,1x,i10/)') + & 'IBO localization (',ltyp,'): IBOs will be stored', + & 'in file locorb.movecs, number ',n1, 'to', n2 + + return + + end + + + +c ================================================================= + + subroutine sym_orth(g_c, g_s, n1, n2) + + implicit none + +#include "global.fh" +#include "mafdecls.fh" +#include "errquit.fh" +#include "stdio.fh" +#include "numerical_constants.fh" + +c Given input vectors collected in matrix C, and a metric S, +c calculate +c C' = C times [transpose(C) S C]**(-1/2) such that +c transpose(C') S C' = unit matrix + +c matrix S is square symmetric, dimension n1, n1 +c matrix C has dimension n1, n2 + +c in the calling routine, matrix C may be allocated with more than +c n2 columns, which is the reason we use ga_copy_patch below + + integer g_c ! in/out + integer g_s ! input + + integer n1, n2 ! input + + integer g_sc, g_csc + + character*(8) pname + + logical debug + +c ----------------------------------------------------------------- + + pname = 'sym_orth' + + debug = .false. .and. ga_nodeid().eq.0 !during development + + if (debug) write(luout,*) 'entering '//pname + + if(.not.ga_create(mt_dbl,n1,n2,'g_sc',0,0,g_sc)) + & call errquit(pname//': error creating g_sc',141, GA_ERR) + + if(.not.ga_create(mt_dbl,n2,n2,'g_csc',0,0,g_csc)) + & call errquit(pname//': error creating g_csc',142, GA_ERR) + + call ga_dgemm('n','n',n1, n2, n1, 1.0d0, + & g_s, g_c, 0.0d0, g_sc) + + call ga_dgemm('t','n',n2, n2, n1, 1.0d0, + & g_c, g_sc, 0.0d0, g_csc) + + call ga_matpow(g_csc, -FP_Half, 1d-12) + + call ga_copy_patch('n', + & g_c, 1, n1, 1, n2, + & g_sc, 1, n1, 1, n2) + + call ga_dgemm('n','n',n1, n2, n2, 1.0d0, + & g_sc, g_csc, 0.0d0, g_c) + + + if (.not. ga_destroy(g_sc)) + & call errquit(pname//': error destroy g_sc',141, GA_ERR) + + if (.not. ga_destroy(g_csc)) + & call errquit(pname//': error destroy g_csc',142, GA_ERR) + + if (debug) write(luout,*) 'leaving '//pname + + return + + end + +c ================================================================= + + subroutine wrap_dgesvd(g_a, g_u, g_vt, values) + + implicit none + +#include "errquit.fh" +#include "global.fh" +#include "mafdecls.fh" + + integer g_a, g_u, g_vt + double precision values(*) + +c wrapper routine for call to dgesvd. This is similar to ga_svd_seq +c in util/ga_it2.F but some of the array dimensions had to be +c adjusted to match the usage case in the virtual IBO construction +c where m < n + + integer n, m, type, l_a, k_a, l_u, k_u, l_vt, k_vt, + & l_work, k_work, lwork, info, nsing + + character*(11) pname + +c ----------------------------------------------------------------- + + pname = 'wrap_dgesvd' + + call ga_inquire(g_a, type, m, n) + + if (type.ne.MT_DBL) then + call errquit (pname//': array type not MT_DBL', + & type, UNKNOWN_ERR) + end if + + nsing = min(m,n) + + if (ga_nodeid() .eq. 0) then + + lwork = 10*max(m,n) + if (.not. ma_push_get(type, m*n, 'svd', l_a, k_a)) + & call errquit(pname//': a',m*n, MA_ERR) + if (.not. ma_push_get(type, m*nsing, 'svd', l_u, k_u)) + & call errquit(pname//': u',m*nsing, MA_ERR) + if (.not. ma_push_get(type, nsing*n, 'svd', l_vt, k_vt)) + & call errquit(pname//': vt',nsing*n, MA_ERR) + if (.not. ma_push_get(type, lwork, 'svd', l_work, k_work)) + & call errquit(pname//': work',lwork, MA_ERR) + + call ga_get(g_a, 1, m, 1, n, dbl_mb(k_a), m) + + call dgesvd('s','s',m,n,dbl_mb(k_a),m,values, + & dbl_mb(k_u),m,dbl_mb(k_vt),nsing, + & dbl_mb(k_work),lwork,info) + + if (info.ne.0) call errquit(pname//':svd failed',info,MEM_ERR) + + call ga_put(g_u, 1, m, 1, nsing, dbl_mb(k_u), m) + call ga_put(g_vt, 1, nsing, 1, n, dbl_mb(k_vt), nsing) + + if (.not. ma_chop_stack(l_a)) call errquit(pname//': svd ma',0, + & MA_ERR) + + end if + + call ga_sync() + call ga_brdcst(1,values,nsing*8,0) + call ga_sync() + + end diff --git a/src/property/localization_driver.F b/src/property/localization_driver.F index b8ff16c6ca..cfe09c94df 100644 --- a/src/property/localization_driver.F +++ b/src/property/localization_driver.F @@ -34,8 +34,8 @@ c ================================================================== #include "apiP.fh" #include "prop.fh" #include "bgj.fh" -#include "util.fh" - +#include "util.fh" +#include "nwc_const.fh" c subroutine arguments: integer rtdb ! [input] rtdb handle @@ -46,17 +46,18 @@ c subroutine arguments: c local GA handles: integer g_uc(4), g_smat, g_sc, g_tran integer g_movecs(2), g_cmo(2), g_temp, g_tmp1, g_tmp2 - integer g_s2, g_s12, g_p12, g_p21 - integer g_ctilde, g_iao, g_mo c MA variables: integer l_c, k_c, l_sc, k_sc, l_eval, k_eval, l_occ, k_occ integer l_dip(3), k_dip(3) + integer l_pop, k_pop, l_list, k_list c other local variables: - integer nclosed(2), nopen(2), nvirt(2), ncore, - & ndens, nbf, nmo, nloc, i + integer loc_opt + + integer nclosed(2), nopen(2), nvirt(2), ncore, nocc, nvir, + & ndens, nbf, nmo, nloc, i, natoms integer info, dbg, ispin @@ -92,13 +93,8 @@ c other local variables: logical int_normalize external int_normalize -c variables for IAO basis - integer - & nbf_iao, nbf_iao_mxprim, nbf_iao_mxang, nbf_iao_mxcont, - & nbf_iao_mxnbf_cn, nbf_iao_mxnbf_ce, nshells_iao, - & minbas - character*(19) pname + character*(3) ltyp c ================================================================== @@ -124,20 +120,52 @@ c ---------------------- write (luout,9000) ! general header if (loctype.eq.1) then call util_print_centered(luout, - & 'Pipek-Mezey localization', 32, .true.) + & 'Pipek-Mezey localization', 25, .true.) elseif (loctype.eq.2) then call util_print_centered(luout, & 'Boys localization', 25, .true.) elseif (loctype.eq.3) then call util_print_centered(luout, - & 'IAO-IBO localization', 20, .true.) + & 'IAO-IBO localization', 25, .true.) else write (luout,*) 'loctype =',loctype write (luout,*) 'LocType value not supported. Exiting.' return end if end if - + +c retrieve localization option from rtdb if it exists + + if (.not. rtdb_get(rtdb, 'prop:loc_opt', mt_int, 1, + & loc_opt )) + & loc_opt = 0 + +c option is currently only implemented for IBOs: + + if (loc_opt.gt.0 .and. (loctype.ne.3)) then + if (oprint) write(luout,'(/1x,a/1x,a)') + & 'localization option not implemented for chosen loc. type', + & 'setting option to 0 ...' + loc_opt = 0 + end if + + if (oprint) then + if (loc_opt .eq.0) then + call util_print_centered(luout, + & '(occupied orbitals)', 25, .false.) + else if (loc_opt.eq.1) then + call util_print_centered(luout, + & '(virtual orbitals)', 25, .false.) + else if (loc_opt.eq.2) then + call util_print_centered(luout, + & '(occ and virt orbitals)', 25, .false.) + else + call errquit(pname//': localization option out of range', + & loc_opt, INPUT_ERR) + end if + end if + + c ------------------------------------------------------- c assemble some data: @@ -160,8 +188,19 @@ c allocate dbl_mb(leval) = orbital energies: if (debug) & write (luout,*) 'loc driver: nclosed,nopen,nvirt',nclosed(1), - & nopen(1), nvirt(1) - + & nopen(1), nvirt(1) + + natoms = 0 + + if (.not. geom_ncent(geom, natoms)) call errquit + & (pname//': geom',0, GEOM_ERR) + + if (natoms.lt.1) call errquit + & (pname//': could not get natoms from rtdb', 911, UNKNOWN_ERR) + + if (natoms.gt.nw_max_atom) call errquit + & (pname//': nw_max_atom too small ', 911, UNKNOWN_ERR) + c Skip localization if we have an open-shell system c (to be extended eventually, localizing alpha and c beta spin orbitals separately): @@ -190,9 +229,15 @@ c dynamically but for now let's make sure we don't get out of c bounds: if (nmo.gt.maxnloc) call - & errquit('loc_driver: maxnloc too small',0, BASIS_ERR) + & errquit('loc_driver: maxnloc too small',0, BASIS_ERR) + +c for the time being, we set these variables here: + + ispin = 1 + nocc = nclosed(1) + nvir = nvirt(ispin) + - c -------------------------------------------- c setup complete. Now proceed to the localization c -------------------------------------------- @@ -471,335 +516,72 @@ c smat not needed anymore if (.not. ga_destroy(g_smat)) call errquit( & 'loc_driver: error destroying g_smat',0, GA_ERR) -c ================================================================ + +c =============================================================== else if (loctype.eq.3) then -c ==================== -c IAO-IBO localization -c ==================== - -c AO Overlap Matrix S1 -> g_smat: - - 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 auxiliary basis overlap S2 -> g_s2 basis needs to defined in the -c input as "iao basis". we will now create the basis here and then -c calculate the overlap and the mixed ao-iao basis overlap S12 - - if (.not. bas_create(minbas, 'iao basis')) - & call errquit(pname//': cannot create iao bas', 86, BASIS_ERR) - - if (.not. bas_rtdb_load(rtdb, geom, minbas, 'iao basis')) - & call errquit(pname//': iao basis not present', 86, BASIS_ERR) - - if (.not. int_normalize(rtdb,minbas)) - & call errquit(pname//': int_normalize failed iao', 86, INT_ERR) - - if (.not. bas_numbf(minbas, nbf_iao)) - & call errquit(pname//': bas error nbf_iao', 86, BASIS_ERR) - - if (.not. bas_nprim_cn_max(minbas,nbf_iao_mxprim)) - & call errquit(pname//': bas error iao mxprim', 86, BASIS_ERR) - - if (.not. bas_high_angular(minbas,nbf_iao_mxang)) - & call errquit(pname//': bas error iao mxang', 86, BASIS_ERR) - - if (.not. bas_ncontr_cn_max(minbas,nbf_iao_mxcont)) - & call errquit(pname//': bas error iao mxcont', 86, BASIS_ERR) - - if (.not. bas_nbf_cn_max(minbas,nbf_iao_mxnbf_cn)) - & call errquit(pname//': bas error iao mcnbf_cn', 86, BASIS_ERR) - - if (.not. bas_nbf_ce_max(minbas,nbf_iao_mxnbf_ce)) - & call errquit(pname//': bas error iao mxnbf_ce', 86, BASIS_ERR) - - if (.not. bas_numcont(minbas,nshells_iao)) - & call errquit(pname//': bas error iao nshells', 86, BASIS_ERR) - -c create overlap for minbas, and the mixed basis-minbas -c overlap S12. -c Then we calculate the projectors P12 and P21 -c via solving S1 P12 = S12 for P12 etc. - - g_s2 = ga_create_atom_blocked(geom, minbas, 'loc:s2') - call ga_zero(g_s2) - call int_1e_ga(minbas, minbas, g_s2, 'overlap', .false.) - - if (.not. ga_create(MT_DBL, nbf, nbf_iao, 'loc:s12', - & nbf, 0, g_s12)) call errquit(pname//': S12',61, GA_ERR) - - call ga_zero(g_s12) - call int_1e_ga(basis, minbas, g_s12, 'overlap', .false.) - - if (.not.ga_duplicate(g_s12, g_p12, 'g_p12')) call - & errquit(pname//': ga_dup P12',62,GA_ERR) - - call ga_copy(g_s12, g_p12) - -#ifdef SCALAPACK - info = ga_solve(g_smat, g_p12) - if (debug) write(luout,*) - & 'info from ga_solve P12:', info -#else - if (debug) write(luout,*) 'call 1 ma_solve' - call ma_solve(g_smat, g_p12) -#endif - - if (.not. ga_create(MT_DBL, nbf_iao, nbf, 'loc:p21', - & nbf_iao, 0, g_p21)) call errquit(pname//': P21',63, GA_ERR) - - call ga_transpose(g_s12, g_p21) - if (debug) write(luout,*) 'P12 assembled' - -c S12 no longer needed - if (.not. ga_destroy(g_s12)) - & call errquit(pname//': ga_destroy failed g_s12',61, GA_ERR) - -c : scalapack throws a 'ga_to_SL: rows/cols error' so we're using -c ma_solve for the time being -c#ifdef SCALAPACK -c info = ga_solve(g_s2, g_p21) -c if (debug) write(luout,*) -c & 'info from ga_solve P21:', info -c#else - if (debug) write(luout,*) 'call 2 ma_solve' - call ma_solve(g_s2, g_p21) -c#endif - -c g_s2 no longer needed - if (.not. ga_destroy(g_s2)) call errquit( - & pname//': error destroying g_s2',0, GA_ERR) - - if (dbg>1) then - call ga_print(g_p12) - call ga_print(g_p21) - end if - -c store P12 * P21 in array g_temp; we'll need it soon - - if (.not.ga_duplicate(g_smat, g_temp, 'g_temp')) call - & errquit(pname//': ga_dup temp', 71, GA_ERR) - - call ga_dgemm('n', 'n', nbf, nbf, nbf_iao, - & 1.0d0, g_p12, g_p21, 0.0d0, g_temp) - - if (debug) write (luout,*) 'p12 p21 -> temp' - -c P21 is no longer needed - - if (.not. ga_destroy(g_p21)) - & call errquit(pname//': ga_destroy failed g_p21',61, GA_ERR) - -c we will localize the doubly occupied orbitals, assuming a -c closed-shell calculation and no designated core -c corbitals. Eventually, this will need to be extended to a loop -c over spins, for unrestricted SCF, and we may want to add -c localization of the virtual space, too - - if (.not. geom_num_core(rtdb, geom, 'ddscf', ncore)) ncore = 0 - - if (debug) write (luout,*) 'ncore = ',ncore - - ispin = 1 - - do i = 1, nclosed(ispin) - iloc(i) = i - end do - nloc = nclosed(ispin) - -c copy the relevant CMOs to g_mo - - if (.not. ga_create(MT_DBL, nbf, nloc, 'loc:g_mo', - & nbf, 0, g_mo)) call errquit(pname//': g_mo',0, GA_ERR) - - call ga_copy_patch('n', - & g_movecs(ispin), 1, nbf, 1, nloc, - & g_mo, 1, nbf, 1, nloc) - - if (debug) write (luout,*) 'movecs -> mo' - -c create C-tilde from Appendix C of Knizia's IBO paper. g_temp -c holds P12 * P21; we won't need it after the next matrix -c multiplication - - if (.not. ga_create(MT_DBL, nbf, nloc, 'loc:ctilde', - & nbf, 0, g_ctilde)) call errquit('loc_driver: sc',0, GA_ERR) - - call ga_dgemm('n', 'n', nbf, nloc, nbf, - & 1.0d0, g_temp, g_mo, 0.0d0, g_ctilde) - - if (.not. ga_destroy(g_temp)) - & call errquit(pname//': ga_destroy failed g_temp',71, GA_ERR) - - if (debug) write (luout,*) 'pre sym_orth 1' - call sym_orth(g_ctilde, g_smat, nbf, nloc) - - if (dbg>1) call ga_print(g_ctilde) - -c now assemble the matrix products to form matrix A from Appendix -c C of the IBO paper. The result will be accumulated in array -c g_iao which are the IAOs in basis 1. The shape of the array is -c the same as P12 - - if (.not.ga_duplicate(g_smat, g_tmp1, 'g_tmp1')) call - & errquit(pname//': ga_dup tmp1', 72, GA_ERR) - - if (.not.ga_duplicate(g_smat, g_tmp2, 'g_tmp2')) call - & errquit(pname//': ga_dup tmp2', 72, GA_ERR) - - if (.not.ga_duplicate(g_smat, g_temp, 'g_temp')) call - & errquit(pname//': ga_dup temp', 72, GA_ERR) - - if (.not.ga_duplicate(g_p12, g_iao, 'g_iao')) call - & errquit(pname//': ga_dup iao', 72, GA_ERR) - - -c C * C(trans) - call ga_dgemm('n', 't', nbf, nbf, nloc, - & 1.0d0, g_mo, g_mo, 0.0d0, g_temp) - -c C * C(trans) * S1 -> tmp1 - call ga_dgemm('n', 'n', nbf, nbf, nbf, - & 1.0d0, g_temp, g_smat, 0.0d0, g_tmp1) - -c tilde-C * tilde-C(trans) - call ga_dgemm('n', 't', nbf, nbf, nloc, - & 1.0d0, g_ctilde, g_ctilde, 0.0d0, g_temp) - -c tilde-C * tilde-C(trans) * S1 -> tmp2 - call ga_dgemm('n', 'n', nbf, nbf, nbf, - & 1.0d0, g_temp, g_smat, 0.0d0, g_tmp2) - -c C * C(trans) * S1 * tilde-C * tilde-C(trans) * S1 - - call ga_dgemm('n', 'n', nbf, nbf, nbf, - & 1.0d0, g_tmp1, g_tmp2, 0.0d0, g_temp) - -c multiply prev. result by P12 and store in g_iao - - call ga_dgemm('n', 'n', nbf, nbf_iao, nbf, - & 1.0d0, g_temp, g_p12, 0.0d0, g_iao) - -c form C * C(trans) -1 - call ga_add_todiag(g_tmp1, -1.0d0) - -c form tilde-C * tilde-C(trans) -1 - call ga_add_todiag(g_tmp2, -1.0d0) - -c multiply the two prev. results with each other, then multiply -c the result with P12 and add to g_iao - - call ga_dgemm('n', 'n', nbf, nbf, nbf, - & 1.0d0, g_tmp1, g_tmp2, 0.0d0, g_temp) - - call ga_dgemm('n', 'n', nbf, nbf_iao, nbf, - & 1.0d0, g_temp, g_p12, 1.0d0, g_iao) - -c g_ctilde no longer needed - - if (.not. ga_destroy(g_ctilde)) - & call errquit(pname//': ga_destroy failed ctilde',62, GA_ERR) - -c delete temp arrays before calling sym_orth, which allocates its -c own temp memory - - if (.not. ga_destroy(g_tmp1)) - & call errquit(pname//': ga_destroy failed g_tmp1',71, GA_ERR) - - if (.not. ga_destroy(g_tmp2)) - & call errquit(pname//': ga_destroy failed g_tmp2',71, GA_ERR) - - if (.not. ga_destroy(g_temp)) - & call errquit(pname//': ga_destroy failed g_temp',71, GA_ERR) -c last step: orthonormalize g_iao - - if (debug) write (luout,*) 'pre sym_orth 2' - call sym_orth(g_iao, g_smat, nbf, nbf_iao) - - call ga_sync - - if (dbg>1) call ga_print(g_iao) - -c check orthogonality of IAOs w.r.t. metric S1 -c g_p12 is used for scratch - - if (.not. ga_create(MT_DBL, nbf_iao, nbf_iao, 'g_temp', - & nbf_iao, 0, g_temp)) call errquit(pname//': temp',73, GA_ERR) - - call ga_dgemm('n', 'n', nbf, nbf_iao, nbf, - & 1.0d0, g_smat, g_iao, 0.0d0, g_p12) - - call ga_dgemm('t', 'n', nbf_iao, nbf_iao, nbf, - & 1.0d0, g_iao, g_p12, 0.0d0, g_temp) - - if (dbg>1) call ga_print(g_temp) - - call ga_add_todiag(g_temp, -1.0d0) - - call ga_maxelt(g_temp, rtemp) - if (master) then - write(luout,'(/1x,a,1x,f18.12)') - & 'IBO loc: largest element in C(iao, trans) S C(iao) -1: ', - & rtemp - write(luout,'(1x,a/1x,a/)') - & 'Significant deviations from zero may indicate', - & 'elevated numerical noise in the IAO generation' - end if - - if (.not. ga_destroy(g_temp)) - & call errquit(pname//': ga_destroy failed g_temp',73, GA_ERR) - -c g_p12 no longer needed - - if (.not. ga_destroy(g_p12)) - & call errquit(pname//': ga_destroy failed g_p12',61, GA_ERR) - -c allocate some memory used in the localization routine: - - if (.not. ma_push_get(mt_dbl, 2*nbf_iao, 'c', l_c, k_c)) - & call errquit(pname//': ma for c', 0, MA_ERR) - - -c jochen: comment -c IBO localization was modeled on the pre-existing PM -c localization code - call localizeIBO(basis, minbas, - & dbl_mb(k_c), - & nloc, iloc, nbf, nmo, g_mo, g_smat, g_iao, nbf_iao, - & dbl_mb(k_eval+(ispin-1)*nbf), dbl_mb(k_occ+(ispin-1)*nbf)) - -c g_mo now holds the localized MOs in the AO basis - -c assuming occupied MO localization, we copy the LMOs into the -c relevant part of movecs - - call ga_copy_patch('n', - & g_mo, 1, nbf, 1, nloc, - & g_movecs(1), 1, nbf, 1, nloc) - - -c deallocate remaining arrays that are no longer needed - - if (.not. ma_pop_stack(l_c)) call errquit( - & pname//': error MA pop c',0, MA_ERR) - - if (.not. ga_destroy(g_iao)) call errquit( - & pname//': error destroying g_iao',0, GA_ERR) - - if (.not. ga_destroy(g_mo)) call errquit( - & pname//': error destroying g_mo',0, GA_ERR) - -c smat not needed anymore - if (.not. ga_destroy(g_smat)) call errquit( - & pname//': error destroying g_smat',0, GA_ERR) - -c destroy minimal basis (iao basis) - if (.not.bas_destroy(minbas)) - & call errquit(pname//'iao bas_destroy failed',0,BASIS_ERR) - +c =================================================== +c IAO-IBO localization (occupied or virtual orbitals) +c =================================================== + + + ispin = 1 + nocc = nclosed(ispin) + nvir = nvirt(ispin) + +c allocate MA array for pairs of MOs. The actual dimensions +c needed are less or equal to 2*nbf + + if (.not. ma_push_get(mt_dbl, 2*nbf, 'c mo', l_c, k_c)) + & call errquit(pname//': ma for c', 0, MA_ERR) + + if (.not. ma_push_get(mt_dbl, natoms, 'pop', l_pop, k_pop)) + & call errquit(pname//': loc:pop', 0, MA_ERR) + + if (.not. ma_push_get(mt_int, natoms, 'list', l_list, k_list)) + & call errquit(pname//': loc:list', 0, MA_ERR) + + if (loc_opt.eq.0) then + ltyp = 'occ' + call ibo_localization(rtdb, geom, ltyp, basis, + & g_movecs(ispin),nocc,nvir, nmo, nbf, natoms, + & dbl_mb(k_eval+(ispin-1)*nbf), + & dbl_mb(k_occ+(ispin-1)*nbf), + & dbl_mb(k_c), + & dbl_mb(k_pop), int_mb(k_list)) + + else if (loc_opt.eq.1) then + ltyp = 'vir' + call ibo_localization(rtdb, geom, ltyp, basis, + & g_movecs(ispin),nocc,nvir, nmo, nbf, natoms, + & dbl_mb(k_eval+(ispin-1)*nbf), + & dbl_mb(k_occ+(ispin-1)*nbf), + & dbl_mb(k_c), + & dbl_mb(k_pop), int_mb(k_list)) + + else if(loc_opt.eq.2) then + ltyp = 'occ' + call ibo_localization(rtdb, geom, ltyp, basis, + & g_movecs(ispin),nocc,nvir, nmo, nbf, natoms, + & dbl_mb(k_eval+(ispin-1)*nbf), + & dbl_mb(k_occ+(ispin-1)*nbf), + & dbl_mb(k_c), + & dbl_mb(k_pop), int_mb(k_list)) + + ltyp = 'vir' + call ibo_localization(rtdb, geom, ltyp, basis, + & g_movecs(ispin),nocc,nvir, nmo, nbf, natoms, + & dbl_mb(k_eval+(ispin-1)*nbf), + & dbl_mb(k_occ+(ispin-1)*nbf), + & dbl_mb(k_c), + & dbl_mb(k_pop), int_mb(k_list)) + + else + call errquit(pname//': loc_opt out of range',loc_opt, + & UNKNOWN_ERR) + + end if ! loc_opt end if ! loctype @@ -849,84 +631,6 @@ c ================================================================== & 10x,28(1h*),/,10x,'*** ORBITAL LOCALIZATION ***',/, & 10x,28(1h*)/) - 9020 format(/,10x,'Open shell calculations not yet implemented', + 9020 format(/,10x,'Spin-unrestricted calculations not yet implemented', & ' for the Localization module. Exiting') end - -c ================================================================= - - subroutine sym_orth(g_c, g_s, n1, n2) - - implicit none - -#include "global.fh" -#include "mafdecls.fh" -#include "errquit.fh" -#include "stdio.fh" -#include "numerical_constants.fh" - -c Given input vectors collected in matrix C, and a metric S, -c calculate -c C' = C times [transpose(C) S C]**(-1/2) such that -c transpose(C') S C = unit matrix - -c matrix S is square symmetric, dimension n1, n1 -c matrix C has dimension n1, n2 - -c in the calling routine, matrix C may be allocated with more than -c n2 columns, which is the reason we use ga_copy_patch below - - integer g_c ! in/out - integer g_s ! input - - integer n1, n2 ! input - - integer g_sc, g_csc - - character*(8) pname - - logical debug - -c ----------------------------------------------------------------- - - pname = 'sym_orth' - - debug = .false. .and. ga_nodeid().eq.0 !during development - - if (debug) write(luout,*) 'entering '//pname - - if(.not.ga_create(mt_dbl,n1,n2,'g_sc',0,0,g_sc)) - & call errquit(pname//': error creating g_sc',141, GA_ERR) - - if(.not.ga_create(mt_dbl,n2,n2,'g_csc',0,0,g_csc)) - & call errquit(pname//': error creating g_csc',142, GA_ERR) - - call ga_dgemm('n','n',n1, n2, n1, 1.0d0, - & g_s, g_c, 0.0d0, g_sc) - - call ga_dgemm('t','n',n2, n2, n1, 1.0d0, - & g_c, g_sc, 0.0d0, g_csc) - - call ga_matpow(g_csc, -FP_Half, 1d-12) - - call ga_copy_patch('n', - & g_c, 1, n1, 1, n2, - & g_sc, 1, n1, 1, n2) - - call ga_dgemm('n','n',n1, n2, n2, 1.0d0, - & g_sc, g_csc, 0.0d0, g_c) - - - if (.not. ga_destroy(g_sc)) - & call errquit(pname//': error destroy g_sc',141, GA_ERR) - - if (.not. ga_destroy(g_csc)) - & call errquit(pname//': error destroy g_csc',142, GA_ERR) - - if (debug) write(luout,*) 'leaving '//pname - - return - - end - - diff --git a/src/property/prop_input.F b/src/property/prop_input.F index 79f41912b9..1d4c5fe842 100644 --- a/src/property/prop_input.F +++ b/src/property/prop_input.F @@ -68,7 +68,7 @@ c ... jochen: integer aoresponse, nfreq, l_fr, k_fr, & velocity, magneticpert, aor_giao, aor_damp, aor_conv, & aor_orbeta, aor_gprime, aor_quad, aor_analyze, - & loc_boys, loc_pm, loc_ibo + & loc_boys, loc_pm, loc_ibo, loc_opt, itmp double precision gamwidth, convergence c c>>> declarations from NWdft @@ -144,7 +144,8 @@ c ... jochen: loc_boys = 1 loc_pm = 1 loc_ibo = 1 -c + loc_opt = 0 +c c read inputs; store in rtdb c do iline = 1,infinity ! infinite loop over input lines @@ -315,15 +316,29 @@ c ... jochen: localization input if (inp_a(test)) then if ( inp_compare(.false., 'boys', test)) then loc_boys = 0 + itmp = 0 + if (inp_i(itmp)) then ! localization option + loc_opt = itmp + end if elseif ( inp_compare(.false., 'ibo', test)) then loc_ibo = 0 + itmp = 0 + if (inp_i(itmp)) then ! localization option + loc_opt = itmp + end if elseif ( inp_compare(.false., 'pm', test)) then loc_pm = 0 + itmp = 0 + if (inp_i(itmp)) then ! localization option + loc_opt = itmp + end if else ! default is PM loc_pm = 0 + loc_opt = 0 end if else loc_pm = 0 + loc_opt = 0 end if c ... jochen end elseif ( inp_compare(.false., 'stoneanalysis', test)) then @@ -687,7 +702,12 @@ c c if (.not. rtdb_put(rtdb, 'prop:loc_ibo', mt_int, 1, $ loc_ibo )) - $ call errquit('prop_input: rtdb_put failed', 0, RTDB_ERR) + $ call errquit('prop_input: rtdb_put failed', 0, RTDB_ERR) +c + if (.not. rtdb_put(rtdb, 'prop:loc_opt', mt_int, 1, + $ loc_opt )) + $ call errquit('prop_input: rtdb_put failed', 0, RTDB_ERR) +c if (.not. rtdb_put(rtdb, 'prop:spinspin', mt_int, 1, $ spinspin )) $ call errquit('prop_input: rtdb_put failed', 0, RTDB_ERR) diff --git a/src/property/prp.F b/src/property/prp.F index 7dfed962d9..dbcf23ff74 100644 --- a/src/property/prp.F +++ b/src/property/prp.F @@ -34,8 +34,8 @@ c 2 nogiao,nogshift,noston,nospinspin,noaoresponse,noboysloc, 3 nopmloc,noiboloc,noelpiso 4 / 1, 1, 1, 1, 1, 1, 1, 1, 1, - 5 1, 1, 1, 1, 1, 1, 1, 1, - 6 1, 1, 1, 1, 1, 1, 1, 1, 1 / + 5 1, 1, 1, 1, 1, 1, 1, 1, 1, + 6 1, 1, 1, 1, 1, 1, 1, 1 / c out =.false. some =.true. @@ -63,6 +63,7 @@ c status=rtdb_get(rtdb,'prop:loc_ibo',MT_INT,1,noiboloc) status=rtdb_get(rtdb,'prop:esp',MT_INT,1,noelp) !done + status=rtdb_get(rtdb,'prop:espiso',MT_INT,1,noelpiso) !done status=rtdb_get(rtdb,'prop:efield',MT_INT,1,noelf) !done status=rtdb_get(rtdb,'prop:efieldgrad',MT_INT,1,noelfg) !done @@ -89,8 +90,8 @@ c Boys localization if (noboysloc.eq.0) then call localization_driver(rtdb, basis, geom, 2) end if -c -c Boys localization + +c IBO localization if (noiboloc.eq.0) then call localization_driver(rtdb, basis, geom, 3) end if