From 19876612be243a8aac466b40b4fdbc92333a6ac6 Mon Sep 17 00:00:00 2001 From: Edoardo Apra Date: Fri, 2 Aug 2013 18:06:50 +0000 Subject: [PATCH] subroutines to generate WFN/WFX files from Alvaro Vazquez-Mayagoitia https://sites.google.com/site/alvarovazquezmayagoitia/goals/codes/nwchem-notes/generator-of-aim-wavefunction-files-nwchem --- QA/doqmtests | 2 +- QA/doqmtests.mpi | 2 +- QA/doqmtests.tcg5 | 2 +- QA/nqetests | 2 + doc/user/property.tex | 14 + src/property/GNUmakefile | 1 + src/property/prop.F | 5 +- src/property/prop_input.F | 9 + src/property/waimfile.F | 1581 +++++++++++++++++++++++++++++++++++++ 9 files changed, 1614 insertions(+), 4 deletions(-) create mode 100644 src/property/waimfile.F diff --git a/QA/doqmtests b/QA/doqmtests index 7364570920..e129c68f7e 100755 --- a/QA/doqmtests +++ b/QA/doqmtests @@ -10,7 +10,7 @@ endif ./runtests.unix procs $np auh2o autosym dft_he2+ h2mp2 h2o hess_h2o prop_h2o pyqa ./runtests.unix procs $np geom_zmatrix rimp2_ne scf_feco5 small_intchk tagcheck testtab ./runtests.unix procs $np h2o_dk u_sodft cosmo_h2o ch5n_nbo h2s_finite startag -./runtests.unix procs $np cosmo_trichloroethene esp esp_uhf dft_bsse dft_s12gh +./runtests.unix procs $np cosmo_trichloroethene esp esp_uhf dft_bsse dft_s12gh c4h4_wfn ./runtests.unix procs $np dplot dft_meta dft_mpwb1khf dft_m05nh2ch3 prop_uhf_h2o ./runtests.unix procs $np et_zn_dimer vectors_rotate sad_ch3hf diff --git a/QA/doqmtests.mpi b/QA/doqmtests.mpi index f38060311f..4c92310ae0 100755 --- a/QA/doqmtests.mpi +++ b/QA/doqmtests.mpi @@ -28,7 +28,7 @@ endif ./runtests.mpi.unix procs $np geom_zmatrix rimp2_ne rimp2_he scf_feco5 small_intchk tagcheck testtab ./runtests.mpi.unix procs $np h2o_dk u_sodft cosmo_h2o ch5n_nbo h2s_finite startag ./runtests.mpi.unix procs $np cosmo_h2o_dft cosmo_h2o_bq be dft_s12gh -./runtests.mpi.unix procs $np cosmo_trichloroethene esp esp_uhf dft_bsse bsse_dft_trimer +./runtests.mpi.unix procs $np cosmo_trichloroethene esp esp_uhf dft_bsse bsse_dft_trimer c4h4_wfn ./runtests.mpi.unix procs $np cosmo_h2cco2 cosmo_h2cco2mg cosmo_h2cco2mg_ecp ./runtests.mpi.unix procs $np cosmo_h3co cosmo_h3co_ecp cosmo_h2cco2na cosmo_h3co_gp ./runtests.mpi.unix procs $np dplot dft_meta dft_mpwb1khf dft_m05nh2ch3 prop_uhf_h2o diff --git a/QA/doqmtests.tcg5 b/QA/doqmtests.tcg5 index cd66b8fd5a..6821f10b73 100755 --- a/QA/doqmtests.tcg5 +++ b/QA/doqmtests.tcg5 @@ -10,7 +10,7 @@ endif ./runtests.tcg5.unix procs $np auh2o autosym dft_he2+ h2mp2 h2o hess_h2o prop_h2o pyqa ./runtests.tcg5.unix procs $np geom_zmatrix rimp2_ne scf_feco5 small_intchk tagcheck testtab ./runtests.tcg5.unix procs $np h2o_dk u_sodft cosmo_h2o ch5n_nbo h2s_finite startag -./runtests.tcg5.unix procs $np cosmo_trichloroethene esp esp_uhf dft_bsse +./runtests.tcg5.unix procs $np cosmo_trichloroethene esp esp_uhf dft_bsse c4h4_wfn ./runtests.tcg5.unix procs $np dplot dft_meta prop_uhf_h2o dft_s12gh ./runtests.tcg5.unix procs $np et_zn_dimer vectors_rotate sad_ch3hf diff --git a/QA/nqetests b/QA/nqetests index a80394a5d4..6126f29cd2 100755 --- a/QA/nqetests +++ b/QA/nqetests @@ -43,6 +43,8 @@ cd ~/nwchem/QA/tests/neda_nbo nwdebug neda_nbo.nw -queue debug -accnt mp15 -procs 2 -time 1999 cd ~/nwchem/QA/tests/ch5n_nbo nwdebug ch5n_nbo.nw -queue debug -accnt mp15 -procs 2 -time 1999 +cd ~/nwchem/QA/tests/c4h4_wfn +nwdebug c4h4_wfn.nw -queue debug -accnt mp15 -procs 2 -time 1999 #--- small tests that should fail! cd ~/nwchem/QA/tests/oh2 nwdebug oh2.nw -queue debug -accnt mp15 -procs 2 -time 1999 diff --git a/doc/user/property.tex b/doc/user/property.tex index f1af03f9f7..15da8558c8 100644 --- a/doc/user/property.tex +++ b/doc/user/property.tex @@ -52,6 +52,7 @@ Each property can be requested by defining one of the following keywords: HYPERFINE SHIELDING [ number_of_atoms atom_list] SPINSPIN [ number_of_pairs pair_list] + AIMFILE ALL \end{verbatim} @@ -100,3 +101,16 @@ string following the \verb+START+ directive. The input deck may be edited to provide additional options to the NBO calculation, (see the NBO user's manual for details.) +\subsection{AIM file} +\label{sec:AIMfile} + +The keyword {\tt aimfile} creates a \verb+.wfn+ file from a given +eigenvector file (i.e., \verb+.movecs+), geometry and a basis set. +The \verb+.wfn+ file summarizes the single determinant wavefunction in order to be analyzed +through a third-party software. +Among some known analyses are the Bader's approach for atoms in molecules (AIM), the electron +localization function (ELF), the Fermi hole function and energy density analysis. +The follow input directive will create the an AIM file \verb+.wfx+ format. +\begin{verbatim} + set prop:nowfx F +\end{verbatim} diff --git a/src/property/GNUmakefile b/src/property/GNUmakefile index a3106d4c4b..3dd6930651 100644 --- a/src/property/GNUmakefile +++ b/src/property/GNUmakefile @@ -68,6 +68,7 @@ task_raman.o \ raman_input.o \ raman.o \ + waimfile.o \ prop_grid.o # OBJ = aoresponse_giao_rhs.o diff --git a/src/property/prop.F b/src/property/prop.F index 5b05403e31..5521f6020e 100644 --- a/src/property/prop.F +++ b/src/property/prop.F @@ -9,7 +9,7 @@ c * implicit none integer rtdb ! [input] - integer nbofile + integer nbofile,aimfile logical status logical hnd_property external hnd_property @@ -42,6 +42,9 @@ c if (rtdb_get(rtdb,'prop:nbofile',MT_INT,1,nbofile)) then if(nbofile.eq.1) call wnbofile(rtdb) endif + if (rtdb_get(rtdb,'prop:aimfile',MT_INT,1,aimfile)) then + if(aimfile.eq.1) call waimfile(rtdb) + endif c c finish ecce property output module c diff --git a/src/property/prop_input.F b/src/property/prop_input.F index 4d8bd66e05..516a7596fa 100644 --- a/src/property/prop_input.F +++ b/src/property/prop_input.F @@ -25,6 +25,7 @@ c integer hypfile, hypopt integer efgfile, efgopt integer nbofile, nboopt + integer aimfile integer dipole integer quadrupole integer octupole @@ -76,6 +77,7 @@ c c c>>> Default property settings. c + aimfile = 0 nbofile = 0 nboopt = 0 efgfile = 0 @@ -391,6 +393,9 @@ c ========================================== dipole = 0 elseif ( inp_compare(.false., 'nbofile', test)) then nbofile = 1 + elseif ( inp_compare(.false., 'aimfile', test)) then + aimfile = 1 + elseif ( inp_compare(.false., 'all', test)) then c ... jochen: read an option (integer value) to define c what NBO options will appear if (inp_i(nboopt)) then @@ -444,6 +449,7 @@ c what EFG options will appear quadrupole = 0 dipole = 0 nbofile = 1 + aimfile = 0 efgfile = 1 hypfile = 1 gshiftfile= 1 @@ -501,6 +507,9 @@ c if (.not. rtdb_put(rtdb, 'prop:nbofile', mt_int, 1, $ nbofile )) $ call errquit('prop_input: rtdb_put failed', 0, RTDB_ERR) + if (.not. rtdb_put(rtdb, 'prop:aimfile', mt_int, 1, + $ aimfile )) + $ call errquit('prop_input: rtdb_put failed', 0, RTDB_ERR) c ... jochen: also write NBO option key to RTDB: if (.not. rtdb_put(rtdb, 'prop:nboopt', mt_int, 1, $ nboopt )) diff --git a/src/property/waimfile.F b/src/property/waimfile.F new file mode 100644 index 0000000000..b02d7d5d85 --- /dev/null +++ b/src/property/waimfile.F @@ -0,0 +1,1581 @@ + subroutine waimfile(rtdb) + implicit none +#include "errquit.fh" +#include "bas.fh" +#include "geom.fh" +#include "global.fh" +#include "mafdecls.fh" +#include "rtdb.fh" +#include "stdio.fh" +#include "inp.fh" +#include "util.fh" +#include "nwc_const.fh" +#include "geomP.fh" +c +c This subroutine was designed to creat a WFN/WFX file +c within a nwchem movecs file. +c The format used is the same as the requeried for +c codes to do topological analyses as Atoms In Molecules +c theories and Electron Localization Function. +c + integer rtdb + integer geom, basis + character*255 vec_file + character*255 titlel_vec, basis_name + character*32 theory +c + integer numcont ! number of mapped contractions + integer icont ! contraction index + integer type ! type (sp/s/p/d/..) +c integer nprim ! no. of primitives per shell + integer ngeno ! no. of contractions + integer sphcart ! 0/1 for cartesian/shperical + integer iset + logical status +c + integer kprim ! sum of primitives per shell + integer primprim ! sum of the total of primitives used + integer sprimprim ! cartesian primitives after transformation from spherical +c + double precision total_charge +c + character*(nw_max_path_len) wfnfile + character*255 title + character*255 title_vec + character*20 scftype + integer nset + integer nbf + integer nmo(2) + integer g_vecs(2) + integer g_hcore(2) + integer g_dens(2) +c + double precision ekin, virial, epot + double precision etot +c + double precision forderco(21) +c + integer nocc(2) + integer noccs +c + integer l_icent, k_icent, l_ricent, k_ricent + integer l_itype, k_itype, l_ritype, k_ritype + integer ls_ritype, ks_ritype + integer ls_ricent, ks_ricent + integer ls_rexp, ks_rexp + integer ls_rcoef, ks_rcoef + integer l_exp, k_exp, l_rexp, k_rexp + integer l_rprim, k_rprim + integer l_eval(2), k_eval(2) + integer l_vecs, k_vecs + integer l_coef, k_coef, l_rcoef, k_rcoef !basis coefficients + integer l_co, k_co + integer ls_co, ks_co + integer l_tprim, k_tprim + integer l_nprim, k_nprim + integer l_occ(2), k_occ(2) +c + integer idum + integer i,j,k,l,n,m + integer naco +c + integer natoms +c + double precision coord(3,nw_max_atom) + double precision xyz(3) + character*16 tag + character*80 buf + double precision qnuc(nw_max_atom) + character*4 symbol(nw_max_atom) + integer atn + character*16 rhf +c + integer rprims(6) + integer rprimc(6) +c functions per type or orbital + data rprims /1, 3, 5, 7, 9, 11/ !wrong + data rprimc /1, 3, 6, 10, 15, 21/ +c + logical nowfx +c + logical movecs_read_header + external movecs_read_header +c + logical movecs_read + external movecs_read +c + integer ga_create_atom_blocked + external ga_create_atom_blocked +c + logical int_normalize + external int_normalize +c + logical waimfile_printwfn, waimfile_printwfx, waimfile_sph2cart + external waimfile_printwfn , waimfile_printwfx, waimfile_sph2cart +c + integer cartbf(-1:5,1:21) + integer sphebf(-1:5,1:11) +c + + character*4 atmsymbol(109) + data atmsymbol/' H ',' HE',' LI',' BE',' B ',' C ', + * ' N ',' O ',' F ',' NE',' NA',' MG', + * ' AL',' SI',' P ',' S ',' CL',' AR', + * ' K ',' CA',' SC',' TI',' V ',' CR', + * ' MN',' FE',' CO',' NI',' CU',' ZN', + * ' GA',' GE',' AS',' SE',' BR',' KR', + * ' RB',' SR',' Y ',' ZR',' NB',' MO', + * ' TC',' RU',' RH',' PD',' AG',' CD', + * ' IN',' SN',' SB',' TE',' I ',' XE', + * ' CS',' BA',' LA',' CE',' PR',' ND', + * ' PM',' SM',' EU',' GD',' TB',' DY', + * ' HO',' ER',' TM',' YB',' LU',' HF', + * ' TA',' W ',' RE',' OS',' IR',' PT', + * ' AU',' HG',' TL',' PB',' BI',' PO', + * ' AT',' RN',' FR',' RA',' AC',' TH', + * ' PA',' U ',' NP',' PU',' AM',' CM', + * ' BK',' CF',' ES',' FM',' MD',' NO', + * ' LR',' RF',' DB',' SG',' BH',' HS', + * ' MT'/ +c L + cartbf(-1,1) = 1 + cartbf(-1,2) = 2 + cartbf(-1,3) = 3 + cartbf(-1,4) = 4 +c 1 S + cartbf(0, 1) = 1 +c 2 PX + cartbf(1, 1) = 2 +c 3 PY + cartbf(1, 2) = 3 +c 4 PZ + cartbf(1, 3) = 4 +c 5 DXX + cartbf(2, 1) = 5 +c 8 DXY + cartbf(2, 2) = 8 +c 9 DXZ + cartbf(2, 3) = 9 +c 6 DYY + cartbf(2, 4) = 6 +c 10 DYZ + cartbf(2, 5) = 10 +c 7 DZZ + cartbf(2, 6) = 7 +c 11 FXXX + cartbf(3, 1) = 11 +c 14 FXXY + cartbf(3, 2) = 14 +c 15 FXXZ + cartbf(3, 3) = 15 +c 17 FXYY + cartbf(3, 4) = 17 +c 20 FXYZ + cartbf(3, 5) = 20 +c 18 FXZZ + cartbf(3, 6) = 18 +c 12 FYYY + cartbf(3, 7) = 12 +c 16 FYYZ + cartbf(3, 8) = 16 +c 19 FYZZ + cartbf(3, 9) = 19 +c 13 FZZZ + cartbf(3,10) = 13 +c 21 GXXXX + cartbf(4, 1) = 21 +c 24 GXXXY + cartbf(4, 2) = 24 +c 25 GXXXZ + cartbf(4, 3) = 25 +c 30 GXXYY + cartbf(4, 4) = 30 +c 33 GXXYZ + cartbf(4, 5) = 33 +c 31 GXXZZ + cartbf(4, 6) = 31 +c 26 GXYYY + cartbf(4, 7) = 26 +c 34 GXYYZ + cartbf(4, 8) = 34 +c 35 GXYZZ + cartbf(4, 9) = 35 +c 28 GXZZZ + cartbf(4,10) = 28 +c 22 GYYYY + cartbf(4,11) = 22 +c 27 GYYYZ + cartbf(4,12) = 27 +c 32 GYYZZ + cartbf(4,13) = 32 +c 29 GYZZZ + cartbf(4,14) = 29 +c 23 GZZZZ +c cartbf(4,15) = 'zzzz' + cartbf(4,15) = 23 + +c 56 HXXXXX (500) +c cartbf(5, 1) = 'xxxxx' + cartbf(5, 1) = 56 +c 55 HXXXXY (410) +c cartbf(5, 2) = 'xxxxy' + cartbf(5, 2) = 55 +c 54 HXXXXZ (401) +c cartbf(5, 3) = 'xxxxz' + cartbf(5, 3) = 54 +c 53 HXXXYY (320) +c cartbf(5, 4) = 'xxxyy' + cartbf(5, 4) = 53 +c 52 HXXXYZ (311) +c cartbf(5, 5) = 'xxxyz' + cartbf(5, 5) = 52 +c 51 HXXXZZ (302) +c cartbf(5, 6) = 'xxxzz' + cartbf(5, 6) = 51 +c 50 HXXYYY (230) +c cartbf(5, 7) = 'xxyyy' + cartbf(5, 7) = 50 +c 49 HXXYYZ (221) +c cartbf(5, 8) = 'xxyyz' + cartbf(5, 8) = 49 +c 48 HXXYZZ (212) +c cartbf(5, 9) = 'xxyzz' + cartbf(5, 9) = 48 +c 47 HXXZZZ (203) +c cartbf(5,10) = 'xxzzz' + cartbf(5, 10) = 47 +c 46 HXYYYY (140) +c cartbf(5,11) = 'xyyyy' + cartbf(5, 11) = 46 +c 45 HXYYYZ (131) +c cartbf(5,12) = 'xyyyz' + cartbf(5, 12) = 45 +c 44 HXYYZZ (122) +c cartbf(5,13) = 'xyyzz' + cartbf(5, 13) = 44 +c 43 HXYZZZ (113) +c cartbf(5,14) = 'xyzzz' + cartbf(5, 14) = 43 +c 42 HXZZZZ (104) +c cartbf(5,15) = 'xzzzz' + cartbf(5, 15) = 42 +c 41 HYYYYY (050) +c cartbf(5,16) = 'yyyyy' + cartbf(5, 16) = 41 +c 40 HYYYYZ (041) +c cartbf(5,17) = 'yyyyz' + cartbf(5, 17) = 40 +c 39 HYYYZZ (032) +c cartbf(5,18) = 'yyyzz' + cartbf(5, 18) = 39 +c 38 HYYZZZ (023) +c cartbf(5,19) = 'yyzzz' + cartbf(5, 19) = 38 +c 37 HYZZZZ (014) +c cartbf(5,20) = 'yzzzz' + cartbf(5, 20) = 37 +c 36 HZZZZZ (005) +c cartbf(5,21) = 'zzzzz' + cartbf(5, 21) = 36 + + sphebf(-1,1) = 1 + sphebf(-1,2) = 2 + sphebf(-1,3) = 3 + sphebf(-1,4) = 4 + sphebf(0, 1) = 1 + sphebf(1, 1) = 2 + sphebf(1, 2) = 3 + sphebf(1, 3) = 4 + sphebf(2, 1) = 5 + sphebf(2, 2) = 6 + sphebf(2, 3) = 7 + sphebf(2, 4) = 8 + sphebf(2, 5) = 9 + sphebf(3, 1) = 10 + sphebf(3, 2) = 11 + sphebf(3, 3) = 12 + sphebf(3, 4) = 13 + sphebf(3, 5) = 14 + sphebf(3, 6) = 15 + sphebf(3, 7) = 16 + sphebf(4, 1) = 17 + sphebf(4, 2) = 18 + sphebf(4, 3) = 19 + sphebf(4, 4) = 20 + sphebf(4, 5) = 21 + sphebf(4, 6) = 22 + sphebf(4, 7) = 23 + sphebf(4, 8) = 24 + sphebf(4, 9) = 25 + sphebf(5, 1) = 26 + sphebf(5, 2) = 37 + sphebf(5, 3) = 38 + sphebf(5, 4) = 39 + sphebf(5, 5) = 40 + sphebf(5, 6) = 41 + sphebf(5, 7) = 42 + sphebf(5, 8) = 43 + sphebf(5, 9) = 44 + sphebf(5, 10) = 45 + sphebf(5, 11) = 46 + +c> open wfn file +c + buf = ' ' + write(buf,*) ' WFX file creator ' + write(6,*) + write(6,*) + call util_print_centered(6,buf,40,.true.) + write(6,*) +c +c> get title + if(.not. rtdb_cget(rtdb, 'title', 1, title)) title = 'NWChem Job ' + +c> wfn / wfx + if (.not. rtdb_get(rtdb, 'prop:nowfx',mt_log, 1, nowfx)) + $ nowfx=.true. +c + +c> load geometry and symmetry info + if (.not. geom_create(geom, 'geometry')) + $ call errquit('waimfile : geom_create?', 0, GEOM_ERR) + if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) + $ call errquit('waimfile : no geometry ', 0, RTDB_ERR) + if (.not.geom_ncent(geom, natoms)) + $ call errquit('waimfile : error code = ',2,GEOM_ERR) + if (.not. geom_nuc_charge(geom,total_charge)) + $ call errquit('waimfile: failed to read total charge', + & 0, RTDB_ERR) + +c> geometry + do i=1,natoms + if (.not.geom_cent_get(geom, i, tag, + $ xyz, qnuc(i))) + $ call errquit('waimfile : error code = ',3,GEOM_ERR) + if (tag(1:2) .eq. 'bq' .or. + $ ((tag(1:1).eq.'x') .and. (tag(2:2).ne.'e'))) then ! X but not Xe + symbol(i) = ' Bq ' + else + symbol(i) = atmsymbol(int(qnuc(i))) + endif + coord(1,i)=xyz(1) + coord(2,i)=xyz(2) + coord(3,i)=xyz(3) + enddo + + +c> load the basis set and get info about it +c + if (.not. bas_create(basis, 'ao basis')) + $ call errquit('waimfile : bas_create?', 0, BASIS_ERR) + if (.not. bas_rtdb_load(rtdb, geom, basis, 'ao basis')) then + if (.not. bas_rtdb_load(rtdb, geom, basis, 'mo basis')) + $ call errquit('waimfile : no mo or ao basis set', 0, + & RTDB_ERR) + endif +c + if (.not. int_normalize(rtdb,basis)) call errquit + & ('waimfile : int_normalize failed',911, INT_ERR) +c +c> READ vectors from movecs file +c + if (.not. rtdb_cget(rtdb, 'prop:vectors', 1, + $ vec_file)) then + call util_file_name('movecs', .false.,.false.,vec_file ) + endif + call util_file_name_resolve(vec_file, .false.) +c> movecs + if (.not. movecs_read_header(vec_file, title_vec, basis_name, + & scftype, nbf, nset, nmo, 2)) + & call errquit('waimfile : basis set error:', 86, basis_err) +c +c + if (.not. ma_push_get(mt_dbl,nbf,'occ',l_occ(1),k_occ(1))) + & call errquit('waimfile : failed to allocate orb.',nbf,MA_ERR) + if (.not. ma_push_get(mt_dbl,nbf,'eigenval',l_eval(1),k_eval(1))) + & call errquit('waimfile : failed to allocate orb.',nbf,MA_ERR) + if (.not. ma_push_get(mt_dbl,nbf,'movecs_read', + $ l_vecs,k_vecs)) + $ call errquit('waimfile : ma failed', nbf, MA_ERR) + +c if (nset.eq.2) then + if (.not. ma_push_get(mt_dbl,nbf,'occ',l_occ(2),k_occ(2))) + & call errquit('waimfile : failed to allocate orb.',nbf,MA_ERR) + if (.not. ma_push_get(mt_dbl,nbf,'eigenval', + & l_eval(2),k_eval(2))) + & call errquit('waimfile : failed to allocate orb.',nbf,MA_ERR) +c endif + + g_vecs(1) = ga_create_atom_blocked(geom, basis, 'vecs1') + call ga_zero(g_vecs(1)) + if (nset.eq.2) then + g_vecs(2) = ga_create_atom_blocked(geom, basis, 'vecs2') + call ga_zero(g_vecs(2)) + endif +c +c> determine orb occupation + nocc(1) = 0 + nocc(2) = 0 + do iset = 1, nset + if (.not. movecs_read(vec_file, iset, dbl_mb(k_occ(iset)), + $ dbl_mb(k_eval(iset)), g_vecs(iset))) call errquit + $ ('waimfile fragment: failed read fragment MOs',0, + $ INPUT_ERR) + do i=1,nbf +c> count electrons + if ( dbl_mb(k_occ(iset)+i-1).gt.1.0d-1) + $ nocc(iset) = nocc(iset) + 1 + enddo + enddo + +c> total occupation + noccs = nocc(1) + if (nset.gt.1) noccs = noccs + nocc(2) + + +c> total of contractions + if (.not.bas_numcont(basis,numcont)) + $ call errquit('waimfile : error code = ',2,BASIS_ERR) + +c> expand basis representation + + if (.not. ma_push_get(mt_int,numcont,'nprim',l_nprim,k_nprim)) + $ call errquit('waimfile : ma failedi nprim', numcont, MA_ERR) + if (.not. ma_push_get(mt_int,numcont,'tprim', l_tprim, k_tprim)) + $ call errquit('waimfile : ma failed tprim', numcont, MA_ERR) + if (.not. ma_push_get(mt_int,numcont,'rprim',l_rprim,k_rprim)) + $ call errquit('waimfile : ma failed rprim', numcont, MA_ERR) + if (.not. ma_push_get(mt_int,numcont,'icent',l_icent,k_icent)) + & call errquit('waimfile : failed icent icent.',numcont,MA_ERR) +c + primprim = 0 + kprim = 0 + do icont = 1, numcont + if (.not. bas_continfo(basis, icont, type, + $ int_mb(k_nprim+icont-1), ngeno,sphcart)) + $ call errquit('waimfile : error code = ',3,BASIS_ERR) +c + if (sphcart.eq.1) then +c write (6,*) 'SPHERICAL BASIS WARNING' + int_mb(k_rprim+icont-1) = rprims(type+1) + else + int_mb(k_rprim+icont-1) = rprimc(type+1) + endif +c + int_mb(k_tprim+icont-1) = + $ int_mb(k_rprim+icont-1)*int_mb(k_nprim+icont-1) + primprim = primprim + int_mb(k_tprim+icont-1) +c + if (.not. bas_cn2ce(basis, icont, i)) + $ call errquit('waimfile : error code = ',5,BASIS_ERR) +c + int_mb(k_icent+icont-1) = i + kprim = kprim + int_mb(k_nprim+icont-1) + enddo +c +c* do icont = 1, numcont +c* write(6,*) 'debug icont ', icont +c* write(6,*) 'debug nprim ', int_mb(k_nprim+icont-1) +c* write(6,*) 'debug rprim ', int_mb(k_rprim+icont-1) +c* write(6,*) 'debug tprim ', int_mb(k_tprim+icont-1) +c* enddo +c + if (.not. ma_push_get(mt_dbl,kprim,'exponents', l_exp, k_exp)) + $ call errquit('waimfile : ma failed exp', kprim, MA_ERR) + if (.not. ma_push_get(mt_dbl,kprim,'coeffts', l_coef, k_coef)) + $ call errquit('waimfile : ma failed exp', kprim, MA_ERR) + + kprim = 0 + do icont = 1, numcont + if (.not.bas_get_exponent(basis, icont, + $ dbl_mb(k_exp + kprim))) + $ call errquit('waimfile : failed bas_get_exponent',0, + $ BASIS_ERR) + + if (.not.bas_get_coeff(basis, icont, dbl_mb(k_coef + kprim))) + $ call errquit('waimfile : failed bas_get_coeff',0, + & BASIS_ERR) + + kprim = kprim + int_mb(k_nprim+icont-1) + enddo +c + if (.not. ma_push_get( mt_dbl, primprim, 'rexponents', + $ l_rexp, k_rexp)) + $ call errquit('waimfile : ma failed exp', primprim, MA_ERR) + if (.not. ma_push_get( mt_dbl, primprim, 'rcoefficients', + $ l_rcoef, k_rcoef)) + $ call errquit('waimfile : ma failed exp', primprim, MA_ERR) + if (.not. ma_push_get( mt_int, primprim, 'ricent', + $ l_ricent,k_ricent)) + & call errquit('waimfile : failed icent.',primprim,MA_ERR) + if (.not. ma_push_get( mt_int, primprim, 'ritype', + $ l_ritype,k_ritype)) + & call errquit('waimfile : failed itype.',primprim,MA_ERR) +c +c +c + k = 0 + l = 1 + m = 0 + do icont = 1, numcont +c> determine number type + if (.not. bas_continfo(basis, icont, type, + $ idum, ngeno,sphcart)) + $ call errquit('waimfile : error code = ',3,BASIS_ERR) + n=1 + do i = 1 + k, int_mb(k_rprim+icont-1) + k + do j = 1 , int_mb(k_nprim+icont-1) +c> itype 1:numcount to ritype 1:primprim + if (sphcart.eq.1) then + int_mb(k_ritype+l-1) = sphebf(type,n) + else + int_mb(k_ritype+l-1) = cartbf(type,n) + endif +c> exp 1:kprim to rexp 1:primprim + dbl_mb(k_rexp+l-1) = dbl_mb(k_exp+m+j-1) +c> coeff 1:kprim to rcoeff 1:primprim + dbl_mb(k_rcoef+l-1) = dbl_mb(k_coef+m+j-1) +c> icenter 1:nprim to ricenter 1:primprim + int_mb(k_ricent+l-1) = int_mb(k_icent+icont-1) + l = l + 1 + enddo + n = n + 1 + enddo + k = k + int_mb(k_tprim+icont-1) + m = m + int_mb(k_nprim+icont-1) + enddo !icont +c + if (.not. ma_push_get(mt_dbl,primprim*noccs,'coeffients', + $ l_co, k_co)) + $ call errquit('waimfile : ma failed', primprim*noccs, MA_ERR) + + naco = 0 + do iset = 1, nset + do n = 1, nocc(iset) + call ga_get(g_vecs(iset),1, nbf,n,n,dbl_mb(k_vecs),1) + k = 0 + l = 1 + m = 0 + do icont = 1, numcont + do i = 1 + k, int_mb(k_rprim+icont-1) + k + m = m + 1 + do j=1,int_mb(k_nprim+icont-1) +c> expand eigenvector 1:nbf to 1:primprim + dbl_mb(k_co + primprim*naco + l-1) = + $ dbl_mb(k_vecs+m-1) +c> easy human reading + if ( abs( dbl_mb(k_co + primprim*naco +l-1)).lt.1d-30) + $ dbl_mb(k_co + primprim*naco +l-1) = 0.0 + l = l + 1 + enddo !j + enddo ! i + k = k + int_mb(k_tprim+icont-1) + enddo !icont + naco = naco + 1 + enddo + enddo + +c> WFN files is always in catesian + + if (sphcart.eq.1) then + if (.not. waimfile_sph2cart( primprim, + $ noccs, + $ nocc, + $ nset, + $ k_ritype, + $ k_icent, k_ricent, + $ k_exp, k_rexp, + $ k_coef, k_rcoef, + $ k_co, + $ ks_rexp,ks_rcoef, + $ ks_ricent,ks_ritype, + $ ks_co, + $ sprimprim)) + $ call errquit('waimfile: write error wfn',0, GA_ERR) + endif + + +cccccccccccccccccccccccccc VIRIAL THEOREM VALUES ccccccccccccccccccccc +c +c Recover potential energy (V) by evaluating kinetic energy (T) +c +c E_total = T + V +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + + + if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory)) + $ call errquit('waimfile: theory not specified',555, + & INPUT_ERR) +c + if (.not. rtdb_get(rtdb,theory(1:inp_strlen(theory))//':energy' + $ , MT_DBL,1,etot)) + $ etot = 0.0d0 + + + + + if (etot.ne.0.0d0) then + +c> do kinetic energy + call int_init(rtdb,1,basis) +c + g_dens(1) = ga_create_atom_blocked(geom, basis, 'dens1') +c + call ga_zero(g_dens(1)) + call ga_dgemm('n', 't', nbf, nbf, nocc(1), 1.0d0, g_vecs(1), + $ g_vecs(1), 0.0d0, g_dens(1)) + + g_hcore(1) = ga_create_atom_blocked(geom, basis,'kin') + call ga_zero(g_hcore(1)) + call int_1e_ga(basis, basis, g_hcore(1), 'kinetic', .false.) +c + ekin = ga_ddot(g_dens(1),g_hcore(1)) +c + if (nset.eq.2) then + g_dens(2) = ga_create_atom_blocked(geom, basis, 'dens2') + call ga_zero(g_dens(2)) + call ga_dgemm('n', 't', nbf, nbf, nocc(2), 1.0d0, g_vecs(2), + $ g_vecs(2), 0.0d0, g_dens(2)) + g_hcore(2) = ga_create_atom_blocked(geom, basis,'kin') + call ga_zero(g_hcore(2)) + call int_1e_ga(basis, basis, g_hcore(2), 'kinetic', .false.) +c + ekin = ekin + ga_ddot(g_dens(2),g_hcore(2)) + if (.not. ga_destroy(g_hcore(2))) + $ call errquit('waimfile: destroy g_hcore2',0, GA_ERR) + if (.not. ga_destroy(g_dens(2))) + $ call errquit('waimfile: destroy g_dens2',0, GA_ERR) + else + ekin = 2*ekin + endif + + if (.not. ga_destroy(g_hcore(1))) + $ call errquit('waimfile: destroy g_hcore1',0, GA_ERR) + if (.not. ga_destroy(g_dens(1))) + $ call errquit('waimfile: destroy g_dens1',0, GA_ERR) +c + epot = etot - ekin + virial = abs(epot/ekin) +c + if (ga_nodeid() .eq. 0) then + write (6,*) ' Total Energy E = ', etot + write (6,*) ' Kinetic Energy V = ', ekin + write (6,*) ' Potential Energy T = ', epot + write (6,*) ' Virial ABS(V/T) = ', virial + endif +c + call int_terminate() + else + if (ga_nodeid() .eq. 0) + $ write (6,*) ' Total Energy = Zero (NO FOUND)' + virial = 0.0d0 + endif + + call ga_sync() + +cccccccccccccccccccccccccc VIRIAL THEOREM VALUES end cccccccccccccccccc + +c> wavefunction + + rhf='RHF' + if (nset.gt.1) rhf = 'UHF' + + +ccccccccccccccccccccc Print file cccccccccccccccccccccccc +c> Print .WFN file : waimfile_printwfn +c> Print .WFX file : waimfile_printwfx + + if (ga_nodeid() .eq. 0) then + if (sphcart.eq.1) then + if(nowfx) then + call util_file_name('wfn', .false., .false., wfnfile) + if (.not. waimfile_printwfn(wfnfile,title,natoms,nset,sprimprim, + $ nocc,int_mb(ks_ricent),int_mb(ks_ritype), + $ dbl_mb(ks_co), dbl_mb(ks_co+sprimprim*nocc(1)), + $ dbl_mb(ks_rexp), dbl_mb(ks_rcoef), + $ coord,qnuc,symbol, + $ dbl_mb(k_occ(1)),dbl_mb(k_occ(2)), + $ dbl_mb(k_eval(1)),dbl_mb(k_eval(2)), + $ etot,virial,rhf)) + $ call errquit('waimfile: write error wfn',0, GA_ERR) + else + call util_file_name('wfx', .false., .false., wfnfile) + if (.not. waimfile_printwfx( wfnfile,title, + $ natoms,nset,sprimprim, + $ nocc,int_mb(ks_ricent),int_mb(ks_ritype), + $ dbl_mb(ks_co), dbl_mb(ks_co+sprimprim*nocc(1)), + $ dbl_mb(ks_rexp), dbl_mb(ks_rcoef), + $ coord,qnuc,symbol, + $ dbl_mb(k_occ(1)),dbl_mb(k_occ(2)), + $ dbl_mb(k_eval(1)),dbl_mb(k_eval(2)), + $ etot,virial,rhf,total_charge)) + $ call errquit('waimfile: write error wfn',0, GA_ERR) + endif !nowfx + else + if(nowfx) then + call util_file_name('wfn', .false., .false., wfnfile) + if (.not. waimfile_printwfn( wfnfile,title,natoms,nset,primprim, + $ nocc,int_mb(k_ricent),int_mb(k_ritype), + $ dbl_mb(k_co), dbl_mb(k_co+primprim*nocc(1)), + $ dbl_mb(k_rexp), dbl_mb(k_rcoef), + $ coord,qnuc,symbol, + $ dbl_mb(k_occ(1)),dbl_mb(k_occ(2)), + $ dbl_mb(k_eval(1)),dbl_mb(k_eval(2)), + $ etot,virial,rhf)) + $ call errquit('waimfile: write error wfn',0, GA_ERR) + else + call util_file_name('wfx', .false., .false., wfnfile) + if (.not. waimfile_printwfx( wfnfile,title, + $ natoms,nset,primprim, + $ nocc,int_mb(k_ricent),int_mb(k_ritype), + $ dbl_mb(k_co), dbl_mb(k_co+primprim*nocc(1)), + $ dbl_mb(k_rexp), dbl_mb(k_rcoef), + $ coord,qnuc,symbol, + $ dbl_mb(k_occ(1)),dbl_mb(k_occ(2)), + $ dbl_mb(k_eval(1)),dbl_mb(k_eval(2)), + $ etot,virial,rhf,total_charge)) + $ call errquit('waimfile: write error wfn',0, GA_ERR) + endif !nowfx + endif !sphcart + endif !ga.eq.0 + +ccccccccccccccccccccc Print file END cccccccccccccccccccc +c + + + if (nset.gt.1) then + if (.not. ga_destroy(g_vecs(2))) + $ call errquit('waimfile: destroy g_movecs',0, GA_ERR) + endif + + if (.not. ga_destroy(g_vecs(1))) + $ call errquit('waimfile: destroy g_movecs',0, GA_ERR) + + if(.not.ma_chop_stack(l_occ(1))) + & call errquit('waimfile, ma_chop_stack of l_icent failed',911, + & ma_err) +c +c* if (.not. geom_print(geom)) +c* $ call errquit('property: geom_print ?',0, GEOM_ERR) +c* call util_flush(luout) +c + if (.not.( + & (bas_destroy(basis)) + & .and. + & (geom_destroy(geom)) + & )) + & call errquit + & ('waimfile:error destroying geom and basis handles',911, + & GEOM_ERR) +c + return + end + + logical function waimfile_printwfn(filename,title,natoms, + $ nset,primprim, + $ nocc,ricent,ritype, + $ coa,cob, + $ rexp,rcoef,coord,qnuc,symbol, + $ occa,occb,evala,evalb, + $ etot, virial,rhf) + implicit none + character*80 filename !in + character*80 title !in + integer natoms + integer iset,nset + integer primprim + integer nocc(2) + integer ricent(primprim), ritype(primprim) + integer i,j + integer naco + double precision rexp(primprim),rcoef(primprim) + double precision coord(3,natoms), qnuc(natoms) + double precision coa(primprim,nocc(1)),cob(primprim,nocc(2)) + character*4 symbol(natoms) + double precision occa(*), occb(*) + double precision evala(*), evalb(*) + double precision etot, virial + character*16 rhf +c + integer unitno + parameter (unitno = 38) + + waimfile_printwfn=.false. + write(6,*) ' Name of Final wnf-file ', filename +cc + open(unit=unitno, file=filename, status='unknown', + & form='formatted') + + write(unitno,1001) title(1:80) + write(unitno,1002) nocc(1)+nocc(2),primprim,natoms +ccc write geometry + write(unitno,1003) (symbol(i),i, i, + $ coord(1,i),coord(2,i),coord(3,i), + $ qnuc(i),i=1,natoms) + write(unitno,1004) (ricent(i),i=1,primprim) + write(unitno,1005) (ritype(i),i=1,primprim) + write(unitno,1006) (rexp(i), i=1,primprim) +c +c>>occupation and orb. ener. +c>>coeff X eigenvec + naco=1 + do i = 1, nocc(1) + write(unitno,1007) i,occa(i),evala(i) + write(unitno,1008) (coa(j,naco) *rcoef(j),j=1, primprim) + naco = naco + 1 + enddo + if(nset.eq.2) then + naco=1 + do i = 1, nocc(2) + write(unitno,1007) i,occb(i),evalb(i) + write(unitno,1008) (cob(j,naco) *rcoef(j),j=1, primprim) + naco = naco + 1 + enddo + endif + +c>>end message + write(unitno,1009) +c>> restricted or unretricted banner + write(unitno,1010) rhf,etot, virial +c + close(unit=unitno, status='keep') + + waimfile_printwfn=.true. + +ccc>>>>> WFN Format <<<<<<<<<<<<<<<<<<< + 1001 format (1A80) + 1002 format ('GAUSSIAN',10X,I5,' MOL ORBITALS',1X,I6,' PRIMITIVES', + & 4X,I5,' NUCLEI') + 1003 format (A4,I4,4X,'(CENTRE',I3,')',1X,3F12.8,' CHARGE =',F5.1) + 1004 format ('CENTRE ASSIGNMENTS',2X,20I3) + 1005 format ('TYPE ASSIGNMENTS',4X,20I3) + 1006 format ('EXPONENTS',1X,1P,5E14.7) + 1007 format ('MO',I3,21X,'OCC NO = ',F12.8,' ORB. ENERGY =',F13.8 ) + 1008 format (1P,5E16.8) + 1009 format ('END DATA') + 1010 format (A8,' ENERGY =',F20.10,' VIRIAL(-V/T) =',F13.8) +ccc>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + end + logical function waimfile_printwfx(filename,title,natoms, + $ nset,primprim, + $ nocc,ricent,ritype, + $ coa,cob, + $ rexp,rcoef,coord,qnuc,symbol, + $ occa,occb,evala,evalb, + $ etot, virial,rhf,total_charge) + implicit none + character*80 filename !in + character*80 title !in + integer natoms + integer iset,nset + integer primprim + integer nocc(2) + integer ricent(primprim), ritype(primprim) + integer i,j + integer naco + double precision rexp(primprim),rcoef(primprim) + double precision coord(3,natoms), qnuc(natoms) + double precision coa(primprim,nocc(1)),cob(primprim,nocc(2)) + character*4 symbol(natoms) + character*30 string_var + double precision occa(*), occb(*) + double precision evala(*), evalb(*) + double precision etot, virial + double precision total_charge + character*16 rhf +c + integer unitno + parameter (unitno = 38) + + waimfile_printwfx=.false. + + + write(6,*) ' Name of Final wnf-file ', filename +c + open(unit=unitno, file=filename, status='unknown', + & form='formatted') + + write(unitno,2001) title(1:80) + write(unitno,2004) ' GTO' + write(unitno,2007) natoms + write(unitno,2010) nocc(1)+nocc(2) + write(unitno,2013) 0 + if (nset.eq.1) then + total_charge=total_charge-2*nocc(1) + write(unitno,2016) total_charge + else + total_charge=total_charge-(nocc(1)+nocc(2)) + write(unitno,2016) total_charge + endif + if (nset.eq.1) then + write(unitno,2019) 2*nocc(1) + write(unitno,2022) nocc(1) + write(unitno,2025) nocc(1) + i = 1 + write(unitno,2028) i + else + write(unitno,2019) nocc(1)+nocc(2) + write(unitno,2022) nocc(1) + write(unitno,2025) nocc(2) + write(unitno,2028) (nocc(1)-nocc(2))+1 + endif + write(unitno,2031) 0 + +c write geometry + write(unitno,2034) + do i=1,natoms + write(string_var, '(I30)') i + write(unitno,2035) symbol(i),adjustl(string_var) + enddo + write(unitno,2036) + + write(unitno,2037) + write(unitno,2038) (int(qnuc(i)),i=1,natoms) + write(unitno,2039) + + write(unitno,2040) + write(unitno,2041) (qnuc(i),i=1,natoms) + write(unitno,2042) + + write(unitno,2043) + write(unitno,2044) (coord(1,i),coord(2,i),coord(3,i), + $ i=1,natoms) + + write(unitno,2045) + write(unitno,2046) + write(unitno,2047) primprim + write(unitno,2048) + write(unitno,2049) + write(unitno,2050) (ricent(i),i=1,primprim) + write(unitno,2051) + write(unitno,2052) + write(unitno,2050) (ritype(i),i=1,primprim) + write(unitno,2053) + write(unitno,2054) + write(unitno,2055) (rexp(i), i=1,primprim) + write(unitno,2056) + + +c>>Occupation Numbers + write(unitno,2057) + write(unitno,2041) (occa(i),i=1,nocc(1)) + if(nset.eq.2) write(unitno,2041) (occb(i),i=1,nocc(2)) + write(unitno,2058) +c Orbital Energies + write(unitno,2059) + write(unitno,2041) (evala(i),i=1,nocc(1)) + if(nset.eq.2) write(unitno,2041) (evalb(i),i=1,nocc(2)) + write(unitno,2060) +c Orbital types + write(unitno,2061) + if (nset.eq.1) then + do i = 1,nocc(1) + write(unitno,2062) + enddo + else + do i = 1,nocc(1) + write(unitno,2063) + enddo + do i = 1,nocc(2) + write(unitno,2064) + enddo + endif + write(unitno,2065) + + write(unitno,2066) + naco = 1 + do i = 1, nocc(1) + write(unitno,2067) + write(unitno,2047) i + write(unitno,2068) + write(unitno,2069) (coa(j,naco) *rcoef(j),j=1, primprim) + naco = naco + 1 + enddo + if (nset.eq.2) then + do i = 1, nocc(1) + write(unitno,2067) + write(unitno,2047) i+nocc(1) + write(unitno,2068) + write(unitno,2069) (coa(j,naco) *rcoef(j),j=1, primprim) + naco = naco + 1 + enddo + endif + write(unitno,2070) +cEnergy print + write(unitno,2071) + write(unitno,2041) etot + write(unitno,2072) +cVirial print + write(unitno,2073) + write(unitno,2041) virial + write(unitno,2074) + + if (nset.eq.1) then + write(unitno,2075) 'Restricted SCF' + else + write(unitno,2075) 'Unrestricted SCF' + endif + + close(unit=unitno, status='keep') + + waimfile_printwfx=.true. + + + +c>>>> WFX Format <<<<<<<<<<<<<<<< + + 2001 format ('',/,1A80,/,'',/) + + 2004 format ('',/,1A5,/,'',/) + + 2007 format ('',/,I5,/,'',/) + + 2010 format ('',/,I5,/, + $ '',/) + + 2013 format ('',/,I2,/, + $ ' ',/) + + 2016 format ('',/,F5.2,/,'',/) + + 2019 format ('',/,I3,/,'',/) + + 2022 format ('',/,I5,/, + $ '',/) + + 2025 format ('',/,I5,/, + $ '',/) + + 2028 format ('',/,I5,/, + $ '',/) + + 2031 format ('',/,I5,/, + $ '',/) + + 2034 format ('') + 2035 format (A3,A30) + 2036 format ('',/) + + 2037 format ('') + 2038 format (I3) + 2039 format ('',/) + + 2040 format ('') + 2041 format (E20.12) + 2042 format ('',/) + + 2043 format ('') + 2044 format (1X,3E20.12) + 2045 format ('',/) + + 2046 format ('') + 2047 format (I5) + 2048 format ('',/) + + 2049 format ('') + 2050 format (5I20) + 2051 format ('',/) + + 2052 format ('') + 2053 format ('',/) + + 2054 format ('') + 2055 format (5E20.12) + 2056 format ('',/) + + 2057 format ('') + 2058 format ('',/) + + 2059 format ('') + 2060 format ('',/) + + 2061 format ('') + 2062 format (' Alpha and Beta') + 2063 format (' Alpha') + 2064 format (' Beta') + 2065 format ('',/) + + 2066 format ('') + 2067 format ('') + 2068 format ('',/) + 2069 format (4E20.12) + 2070 format ('',/) + + 2071 format ('') + 2072 format ('',/) + + 2073 format ('') + 2074 format ('',/) + + 2075 format ('',/,1A80,/,'',/) + + end + logical function waimfile_sph2cart( primprim, + $ noccs, + $ nocc, + $ nset, + $ k_ritype, + $ k_icent, k_ricent, + $ k_exp, k_rexp, + $ k_coef, k_rcoef, + $ k_co, + $ ks_rexp,ks_rcoef, + $ ks_ricent,ks_ritype, + $ ks_co, + $ sprimprim) + implicit none +#include "mafdecls.fh" +#include "errquit.fh" + integer primprim + integer sprimprim + integer nset,iset + integer numd,numf,numg,numh + integer noccs, nocc(2) + integer k_icent, k_ricent + integer k_ritype + integer k_exp, k_rexp + integer k_coef, k_rcoef + integer k_co + integer ls_co, ks_co + integer ls_ritype, ks_ritype + integer ls_ricent, ks_ricent + integer ls_rexp, ks_rexp + integer ls_rcoef, ks_rcoef + integer i,j,k,l,n,r,h,m,t + integer np,iorb,forb,nsf + integer naco + double precision sqrt3, sqrt5, sqrt7, DUMMY +c: +c: author A. Guevara 7/2013 +c: + waimfile_sph2cart=.false. + + !sqrt3=1.7320508076D+00 + sqrt3=3.0**(1.0/3.0) + + numd=0 + numf=0 + numg=0 + numh=0 + do i = 0,primprim-1 + if (int_mb(k_ritype+i).eq.5) numd=numd+1 + if (int_mb(k_ritype+i).eq.10) numf=numf+1 + if (int_mb(k_ritype+i).eq.17) numg=numg+1 + if (int_mb(k_ritype+i).eq.26) numh=numh+1 + enddo + sprimprim=primprim+numd+3*numf+6*numg+10*numh + if (.not. ma_push_get(mt_dbl,sprimprim,'rexponents', + $ ls_rexp, ks_rexp)) + $ call errquit('waimfile_sph2cart : ma failed exp', + $ sprimprim, MA_ERR) + if (.not. ma_push_get(mt_dbl,sprimprim,'rcoefficients', + $ ls_rcoef, ks_rcoef)) + $call errquit('waimfile_sph2cart : ma failed exp', + $ sprimprim, MA_ERR) + if (.not. ma_push_get(mt_int,sprimprim,'ricent', + $ ls_ricent,ks_ricent)) + &call errquit('waimfile_sph2cart : failed icent.',sprimprim,MA_ERR) + if (.not. ma_push_get(mt_int,sprimprim,'ritype', + $ ls_ritype,ks_ritype)) + & call errquit('waimfile_sph2cart : failed itype.', + $ sprimprim, MA_ERR) + + + if (.not. ma_push_get(mt_dbl,sprimprim*noccs,'coeffients', + $ ls_co, ks_co)) + $ call errquit('waimfile : ma failed ls coeff', + $ sprimprim*noccs, MA_ERR) + + i=0 + t=0 + do while (i