← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/property/wefgfile.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/property/wefgfile.F
 
1
c*****************************************************************************
 
2
      subroutine wefgfile(rtdb,
 
3
     &                    g_munu_rot,
 
4
     &                    nlst,         
 
5
     &                    Ndir_munu,
 
6
     &                    Natoms_munu,
 
7
     &                    atmnr_munu,
 
8
     &                    Nel)
 
9
c*****************************************************************************
 
10
c
 
11
c>>>  Purpose of routine:
 
12
c
 
13
c>>>  This subroutine reads information from the runtime database
 
14
c>>>  and prints it out to a file, name=nbo_fname, which is used
 
15
c>>>  as input to the standalone version of the natural bond orbital
 
16
c>>>  (NBO) analysis program, gennbo.
 
17
c
 
18
c>>>  Dependencies/Effects:
 
19
c
 
20
c     The rtdb is read only and not modified.
 
21
c     No common blocks are created.
 
22
c     A file, (lfn=25, name=nbo_fname), is created.
 
23
c     All variables except for rtdb are local.
 
24
c
 
25
c>>>  Commentary:
 
26
c
 
27
c>>>  Eventually, it may be desirable to have the NBO analysis executed
 
28
c>>>  as a subroutine of NWchem rather than as a separate stand alone
 
29
c>>>  program.  To do this one will have to make a subroutine similar
 
30
c>>>  to the one here, but instead of writing the data to an input file
 
31
c>>>  it will need to be passed to NBO interface subroutines,  (see the
 
32
c>>>  NBO programmers manual for details).  The current strategy was
 
33
c>>>  chosen to reduce the binary size of NWchem (the size of the NBO 
 
34
c>>>  executable is not negligible), and to take advantage of the fact
 
35
c>>>  that gennbo already exists and has been ported to a number of 
 
36
c>>>  platforms.
 
37
c
 
38
*
 
39
* $Id: wefgfile.F 21401 2011-11-04 18:22:41Z bert $
 
40
*
 
41
      implicit none
 
42
#include "errquit.fh"
 
43
#include "bas.fh"
 
44
#include "geom.fh"
 
45
#include "global.fh"
 
46
#include "mafdecls.fh"
 
47
#include "rtdb.fh"
 
48
#include "stdio.fh"
 
49
#include "zora.fh" ! FA-04-24-10
 
50
c Using g_Ci defined in zora.fh, computed in dft_zora_scale() from dft_zora_utils.F
 
51
      integer rtdb  
 
52
c
 
53
c>>>  handles
 
54
c
 
55
 
 
56
      integer lbasis, lchg, lcoef,lcoord, leval, lexp, lgeom, lind1, 
 
57
     $    lind2, locc, lprim, lptr, lscr, ltags, ltype,
 
58
     &    lncomp
 
59
      integer g_bo(2), g_fock(2), g_movecs(2),  g_over, g_scr
 
60
c
 
61
c>>>  indices
 
62
c
 
63
      integer ichg, icoef, icoord, ieval, iexp, iind1, iind2, iocc,
 
64
     $    iprim, iptr, iscr, itags, itype,
 
65
     &    incomp
 
66
c
 
67
c>>>  basis set objects
 
68
c
 
69
      integer atn, charge, maxL, nat, nbf, ncont, ngeno, nmo, 
 
70
     $     nprim, sphcart, type
 
71
      character*2   sym, elem
 
72
      character*255 name, title, trans
 
73
c
 
74
c>>>  things related to the NBO input file
 
75
c
 
76
      integer lfn
 
77
      character*255 nbo_fname
 
78
c
 
79
c>>>  MO vector and calculation-type objects
 
80
c
 
81
      character*4   spintype
 
82
      character*255 bas_vecs, movecs_in, title_vecs
 
83
      character*20 scftype_vecs
 
84
      integer  lenocc, nbf_vecs, nset, nmo_vecs(2)
 
85
 
 
86
c
 
87
c>>>  Functions
 
88
c
 
89
      logical int_normalize
 
90
      logical  movecs_read, restricted, movecs_read_header
 
91
      integer  inp_strlen, ga_create_atom_blocked
 
92
 
 
93
      external int_normalize
 
94
      external inp_strlen, ga_create_atom_blocked, movecs_read, 
 
95
     $     movecs_read_header
 
96
 
 
97
c
 
98
c>>>  miscellaneous variables
 
99
c
 
100
      character*5 cname(5)
 
101
      integer i,iat,iat1,icol,icont, 
 
102
     &        ilo,ihi,index,ix,iy,iz,j,
 
103
     $        l, len
 
104
      logical sshell, pshell, dshell, fshell, gshell
 
105
      data cname / 'CS = ', 'CP = ', 'CD = ', 'CF = ', 'CG = ' /
 
106
      integer efgopt
 
107
      integer l_Ci,k_Ci, ! FA-scaling-zora
 
108
     &        l_qscal,k_qscal,noc(2),indq,count_qscal,indx_inc
 
109
      integer ipol,noc1,nocc,nelec
 
110
      logical lzora
 
111
c     Assuming noc1=noc2 not true if (total spin .ne. 0)
 
112
c              (noc_alpha = noc_beta)
 
113
      character*11 fname_efg
 
114
      integer nocc_a,nocc_b,idir,jlo,jhi
 
115
      integer l_buf,k_buf,indx      
 
116
      integer g_munu_rot,nlst,
 
117
     &        Ndir_munu,Natoms_munu,Nel,
 
118
     &        atmnr_munu(Natoms_munu)
 
119
      double precision threshold4data
 
120
      data threshold4data /1e-40/ ! Threshhold for output data
 
121
      if (ga_nodeid().eq.0)
 
122
     & write(*,*) 'Using data threshold :',threshold4data
 
123
c     --> To store ONLY munu principal components xx,yy,zz
 
124
c    Structure of g_munu_rot:
 
125
c    n=nbf  Nr. of basis functions
 
126
c    val_11_1_1 val_22_1_1 ... val_nn_1_1 val_21_1_1 
 
127
c                                         val_31_1_1 val_32_1_1
 
128
c                                         val_41_1_1 val_42_1_1 val_43_1_1 
 
129
c                                         ...
 
130
c    val_11_2_1 val_22_2_1 ... val_nn_2_1 val_21_2_1 
 
131
c                                         val_31_2_1 val_32_2_1
 
132
c                                         val_41_2_1 val_42_2_1 val_43_2_1 
 
133
c                                         ...
 
134
c    ...
 
135
c    val_11_1_2 val_22_1_2 ... val_nn_1_2 val_21_1_2 
 
136
c                                         val_31_1_2 val_32_1_2
 
137
c                                         val_41_1_2 val_42_1_2 val_43_1_2
 
138
c                                         ...
 
139
c    val_11_2_2 val_22_2_2 ... val_nn_2_2 val_21_2_2 
 
140
c                                         val_31_2_2 val_32_2_2
 
141
c                                         val_41_2_2 val_42_2_2 val_43_2_2
 
142
c                                         ...
 
143
c    Notation.- val_PQ_R_S
 
144
c               PQ, unique munu indices
 
145
c               R, direction index = 1,2,3 (in this order) 
 
146
c               S, atom number  R=1,nlist
 
147
c               Correspondence with EFG eigenvalues:
 
148
c               V33 V22 V11
 
149
c               -> If multiplied with corresponding 
 
150
c                  scaling density matrix it produces
 
151
c                  the total EFG eigenvalues in the same order
 
152
c                  as NWChem standard output.            
 
153
c               nlist total number of atoms selected in
 
154
c               keyword of input script: efieldgradZ4 [nlist] ...
 
155
c                    
 
156
 
 
157
c ======= reading from rtdb (nocc_a,nocc_b) === START
 
158
         if (.not. rtdb_get(rtdb, 'prop:Nocc_a',mt_int,
 
159
     $                      1,nocc_a))
 
160
     $      call errquit('prop_input-EFGZ4-Nocc_a: rtdb_put failed',
 
161
     $                   555, RTDB_ERR)
 
162
         if (.not. rtdb_get(rtdb, 'prop:Nocc_b',mt_int,
 
163
     $                      1,nocc_b))
 
164
     $      call errquit('prop_input-EFGZ4-Nocc_b: rtdb_put failed',
 
165
     $                   555, RTDB_ERR)
 
166
c ======= reading from rtdb (nocc_a,nocc_b) === END
 
167
c ------ Lines below are not necessary ------ START
 
168
       lfn = 25
 
169
c      if(ga_nodeid().eq.0) then
 
170
c        call util_file_name('hyp', .false., .false., nbo_fname)
 
171
c        open(lfn,file=nbo_fname,status='unknown')
 
172
c        len = inp_strlen(nbo_fname)
 
173
c        write(LUout,9250) nbo_fname(1:len)
 
174
c      endif 
 
175
c ------ Lines below are not necessary ------ END
 
176
 
 
177
c ... Fredy: read EFG option key fron RTDB.
 
178
      efgopt = 0
 
179
      if (.not.rtdb_get(rtdb, 'prop:efgopt', mt_int, 1,efgopt)) 
 
180
     &    call errquit('wefgfile: efgopt RTDB failed',0, RTDB_ERR)
 
181
c ... Fredy end
 
182
c
 
183
c>>>  load geometry and symmetry info
 
184
c
 
185
      if (.not. geom_create(lgeom, 'geometry'))
 
186
     $     call errquit('wnbofile: geom_create?', 0, GEOM_ERR)
 
187
      if (.not. geom_rtdb_load(rtdb, lgeom, 'geometry'))
 
188
     $     call errquit('wnbofile: no geometry ', 0, RTDB_ERR)
 
189
c     
 
190
c>>>  load the basis set and get info about it
 
191
c
 
192
      if (.not. bas_create(lbasis, 'mo basis'))
 
193
     $     call errquit('prop: bas_create?', 0, BASIS_ERR)
 
194
      if (.not. bas_rtdb_load(rtdb, lgeom, lbasis, 'ao basis')) then
 
195
         if (.not. bas_rtdb_load(rtdb, lgeom, lbasis, 'mo basis'))
 
196
     $        call errquit('wnbofile: no mo or ao basis set', 0,
 
197
     &       RTDB_ERR)
 
198
      end if
 
199
      if (.not. bas_name(lbasis, name, trans))
 
200
     $     call errquit('wnbofile: bas_name?', 0, BASIS_ERR)
 
201
c
 
202
c>>>  Create ga_array for movecs.
 
203
c
 
204
      if ( .not. bas_numbf(lbasis,nbf) )
 
205
     $     call errquit('wnbofile: bas_numbf failed', 1, BASIS_ERR)
 
206
c
 
207
c>>>  For now simply set nmo = nbf.
 
208
c
 
209
      nmo = nbf
 
210
c     
 
211
c>>>    Get information on molecular orbital vectors.
 
212
c       
 
213
      call util_file_name('movecs', .false., .false., movecs_in)
 
214
      if(.not.movecs_read_header(movecs_in, title_vecs, bas_vecs,
 
215
     $     scftype_vecs, nbf_vecs, nset, nmo_vecs, 2))
 
216
     $     call errquit('wnbofile:  failed on movecs_read_header',0,
 
217
     &       DISK_ERR)
 
218
 
 
219
      if(nset.eq.1)then
 
220
        restricted=.true.
 
221
        spintype = '    '
 
222
      else 
 
223
        restricted=.false.
 
224
        spintype = 'OPEN' 
 
225
      endif 
 
226
c
 
227
c>>>  Extract title, geometry, nuclear charges, and label info.
 
228
      if (.not. geom_ncent(lgeom, nat))
 
229
     $     call errquit('wnbofile: geom_ncent failed',1, GEOM_ERR)
 
230
c       
 
231
c>>>  Begin $GENNBO and $NBO keylist data.
 
232
c ++++++++++++++++++++++++++++++++++++++++++
 
233
c ++++++++++++ EFG printout ++++++++++ START
 
234
       if (nset.eq.1) then
 
235
        nocc=nocc_a
 
236
        noc(1)=nocc_a
 
237
       else
 
238
        nocc = nocc_a+nocc_b
 
239
        noc(1)=nocc_a
 
240
        noc(2)=nocc_b
 
241
       endif
 
242
       lzora=.false.
 
243
       if (do_zora .and. so_term.eq.0 .and. 
 
244
     &     .not.(not_zora_scale)) lzora=.true.
 
245
       if (lzora) then  ! spin-free-zora
 
246
        if (.not.ma_alloc_get(mt_dbl,nset*nocc,'Ci',l_Ci,k_Ci)) 
 
247
     &     call errquit('wefgfile: ma failed',0,MA_ERR)
 
248
        if (.not.ma_alloc_get(mt_dbl,nocc,'qscal',l_qscal,k_qscal)) 
 
249
     &     call errquit('wefgfile: ma failed',0,MA_ERR)
 
250
 
 
251
        if (ga_nodeid().eq.0)
 
252
     &   write(*,*) '-------g_Ci---------- START'
 
253
        call ga_print(g_Ci)
 
254
        if (ga_nodeid().eq.0)
 
255
     &   write(*,*) '-------g_Ci---------- END'
 
256
 
 
257
        call ga_get(g_Ci,1,nset,1,nocc,dbl_mb(k_Ci),nset)
 
258
        if (nset.eq.1) then
 
259
         indx_inc=1
 
260
        else if (nset.eq.2) then
 
261
         indx_inc=2
 
262
        else
 
263
         write(*,*) "ERROR in wefgfile:: nset=1 or 2 ..."
 
264
         stop
 
265
        endif
 
266
        count_qscal=1
 
267
        do i=1,nset
 
268
         indx=i ! =1 for alpha (ODD nrs) =2 for beta (even nrs)
 
269
         do j=1,noc(i)
 
270
           indq=k_qscal+count_qscal-1
 
271
           dbl_mb(indq)=dbl_mb(k_Ci+indx-1)
 
272
 
 
273
           if (dbl_mb(indq) .ne. 0.0d0)
 
274
     &      dbl_mb(indq)=dbl_mb(indq)-1.0d0
 
275
 
 
276
          if (ga_nodeid().eq.0) then
 
277
           write(*,7) i,j,
 
278
     &                dbl_mb(indq),1.0d0/(1.0d0+dbl_mb(indq)),
 
279
     &                1.0d0/(1.0d0+dbl_mb(indq))**0.5d0
 
280
  7        format('wefgfile::qscal(',i4,',',i4,')=(',f15.8,',',
 
281
     &             f15.8,',',f15.8,')')
 
282
          endif
 
283
          count_qscal=count_qscal+1
 
284
          indx=indx+indx_inc
 
285
         enddo ! end-loop-j
 
286
        enddo ! end-loop-i
 
287
       endif ! end-if-lzora
 
288
 
 
289
      if( ga_nodeid().eq.0) then  ! ------- nodeid0----- START
 
290
c ... now a loop over the EFG components, idir = 1,3.
 
291
c     number of components was given in namelist
 
292
c     ==>  nlst=nbf*(nbf+1)
 
293
        if (.not.ma_alloc_get(mt_dbl,nlst,'buf',
 
294
     &      l_buf,k_buf))
 
295
     &      call errquit('wefgfile-buf: ma failed',0,MA_ERR)
 
296
        do iat1=1,Natoms_munu
 
297
         iat=atmnr_munu(iat1)
 
298
c=====================================
 
299
c ==== Generate filename per atom here
 
300
c ====================================
 
301
         write (fname_efg,1) iat
 
302
 1       format('efgdata.',i0)
 
303
         write(*,*) 'Writing filename=',fname_efg
 
304
         open (lfn, file=fname_efg, status='unknown')
 
305
         write (lfn,
 
306
     &   '(a,i10,a,i10,a,i10,/a,i10,a,i10,a,i10,a,i10,/a,l3,a/a)') 
 
307
     &   '&inputdata nbasis=',nbf, 
 
308
     &   ', nocc='  ,nocc,
 
309
     &   ', nspin=' ,nset,
 
310
     &   ', nocc_a=',nocc_a,
 
311
     &   ', nocc_b=',nocc_b,
 
312
     &   ', nmo='   ,nmo, 
 
313
     &   ', nelec=' ,nocc, 
 
314
     &   ', threshold = 0.0, lzora=',lzora, 
 
315
     &   ', spinorbit= F , label=''EFG'', components=3', 
 
316
     &   ', nbofile=''nbodata'', antisym = F, fullmat = F /'
 
317
          if (lzora) then  ! spin-free-zora
 
318
c     ----move down: g_Ci scaling factors
 
319
           write (lfn,'(a)') 'ZORA MO scaling factors'
 
320
         write (lfn, '(3e25.16)' ) (dbl_mb(k_qscal+i-1), i = 1,nocc)
 
321
          endif
 
322
          do idir=1,3
 
323
           jlo=nlst*Ndir_munu*(iat1-1)+
 
324
     &         nlst*(idir-1)+1
 
325
           jhi=jlo+nlst-1
 
326
           call dcopy(nlst,0.0d0,0,dbl_mb(k_buf),1)
 
327
           call ga_get(g_munu_rot,1,1,jlo,jhi,
 
328
     &                 dbl_mb(k_buf),1)
 
329
            do i=0,nlst-1
 
330
             if (abs(dbl_mb(k_buf+i)) .lt. threshold4data)
 
331
     &         dbl_mb(k_buf+i)=0.0d0
 
332
            enddo
 
333
            write (lfn,'(a,i1)') 'h1mat for direction ',idir
 
334
            write (lfn,'(3e25.16)' ) 
 
335
     &       (dbl_mb(k_buf+i), i = 0,nlst-1)
 
336
          enddo ! end-loop-directions
 
337
          close(lfn) ! close efgdata
 
338
        enddo ! end-loop-atoms
 
339
         if (.not.ma_free_heap(l_buf)) call
 
340
     &      errquit('munu4nbo: ma_free_heap l_buf',0, MA_ERR)
 
341
         if (lzora) then
 
342
          if (.not.ma_free_heap(l_Ci)) call
 
343
     &      errquit('wshldfile: ma_free_heap l_Ci',0, MA_ERR)
 
344
          if (.not.ma_free_heap(l_qscal)) call
 
345
     &    errquit('wgshiftfile: ma_free_heap l_qscal',0, MA_ERR)
 
346
         endif
 
347
       endif ! ------- nodeid0----- START
 
348
       goto 147
 
349
      if (ga_nodeid().eq.0) then ! ------ if-ga_nodeid-1--- START
 
350
        if(.not. rtdb_cget(rtdb, 'title', 1, title)) title = ' '
 
351
        if(.not.Ma_Push_Get(MT_Dbl, nat*3, 'coordinates', 
 
352
     &          lcoord, icoord))
 
353
     $     call errquit('wnbofile: push_get failed',2, MA_ERR)
 
354
        if(.not.Ma_Push_Get(MT_Dbl, nat, 'charges', lchg, ichg))
 
355
     $     call errquit('wnbofile: push_get failed',3, MA_ERR)
 
356
        if(.not.Ma_Push_Get(MT_Byte, nat*16, 'center tags', 
 
357
     &          ltags, itags))
 
358
     $     call errquit('wnbofile: push_get failed',4, MA_ERR)
 
359
        if(.not.geom_cart_get(lgeom, nat, Byte_MB(itags), 
 
360
     &          Dbl_MB(icoord),Dbl_MB(ichg)))
 
361
     $     call errquit('wnbofile: geom_cart_get failed',5, GEOM_ERR)
 
362
        if (restricted) then
 
363
*ga:1:0
 
364
         if (.not. ga_create(MT_DBL, nbf, nmo, 'wnbofile: MOs',
 
365
     $       0, 0, g_movecs(1))) call errquit('wnbofile: MOs', 0,
 
366
     &       GA_ERR)
 
367
         call ga_zero(g_movecs(1))
 
368
        else
 
369
*ga:1:0
 
370
          if (.not. ga_create(MT_DBL, nbf, nmo, 'wnbofile: alpha MOs',
 
371
     $       0, 0, g_movecs(1))) call errquit('wnbofile: alpha MOs', 0,
 
372
     &       GA_ERR)
 
373
          call ga_zero(g_movecs(1))
 
374
*ga:1:0
 
375
          if (.not. ga_create(MT_DBL, nbf, nmo, 'wnbofile: beta MOs',
 
376
     $       0, 0, g_movecs(2))) call errquit('wnbofile: beta MOs', 0,
 
377
     &       GA_ERR)
 
378
          call ga_zero(g_movecs(2))
 
379
        endif
 
380
 
 
381
        open (lfn, file='nbodata.f47_extra', status='unknown')
 
382
c ++++++++++++ EFG printout ++++++++++ END
 
383
c ++++++++++++++++++++++++++++++++++++++++
 
384
        write(lfn,9000) nat, nbf, spintype
 
385
 9000   format('$GENNBO  UPPER  BODM  BOHR  NATOMS=',i5,
 
386
     $         '  NBAS=',i5,2x,a4,'  $END')
 
387
c
 
388
c ... Fredy: here we write more detailed NBO input if requested. 
 
389
        if (efgopt.eq.0) then ! default
 
390
           write(lfn,'(a)')'$NBO     $END'
 
391
        else if (efgopt.eq.1) then ! generate NLMOs, too, and dump info to files
 
392
           write(lfn,'(a/a)')'$NBO BNDIDX NBONLMO=W AONBO=W AONLMO=W ',
 
393
     $          ' NLMOMO=W STERIC FILE=nbodata $END'
 
394
        else if (efgopt.eq.2) then ! generate even more info, incl AO overlap etc.
 
395
           write(lfn,'(a/a)')'$NBO BNDIDX NBONLMO=W AONBO=W AONLMO=W ',
 
396
     $          ' NLMOMO=W AOMO=W SAO=W STERIC FILE=nbodata $END'
 
397
        else ! other options not supported. Fall back to default
 
398
           write(lfn,'(a)')'$NBO     $END'
 
399
        end if
 
400
c ... jochen end
 
401
      endif   ! ------ if-ga_nodeid-1--- END
 
402
c ------- Destroy ga arrays -----------
 
403
c      if (.not. ga_destroy(g_Ci)) call errquit( ! destroy GA - FA
 
404
c     &  'wefgfile: ga_destroy failed ',0, GA_ERR)
 
405
c
 
406
c>>>  End $GENNBO and $NBO keylist data.
 
407
c
 
408
c>>>  Begin $COORD keylist.
 
409
      if( ga_nodeid().eq.0)then
 
410
        write(lfn,9100)title(1:80)
 
411
 9100   format('$COORD',/,a)
 
412
 
 
413
        do iat = 1,nat
 
414
 
 
415
          ix = 3*(iat-1) + icoord 
 
416
          iy = ix+1
 
417
          iz = iy+1
 
418
          index = itags+(iat-1)*16
 
419
 
 
420
          if(.not.geom_tag_to_element(Byte_MB(index),sym,elem,atn))
 
421
cedo     $         call errquit('wnbofile: geom_tag_to_element failed',1)
 
422
     .         then
 
423
          atn=0
 
424
          charge=0
 
425
       endif
 
426
 
 
427
          charge = Dbl_MB(ichg+iat-1)
 
428
          write(lfn,9150)atn, charge, Dbl_MB(ix), Dbl_MB(iy), Dbl_MB(iz)
 
429
 
 
430
        enddo 
 
431
 
 
432
 9150   format(i3,2x,i3,2x,3f20.10)
 
433
        write(lfn,9200)
 
434
      endif
 
435
c
 
436
c>>>  End $COORD keylist.
 
437
c
 
438
c>>>  Begin $BASIS keylist.
 
439
c
 
440
 
 
441
      if(.not.bas_numcont(lbasis,ncont))
 
442
     $     call errquit('wnbofile: bas_numcont failed',0, BASIS_ERR)
 
443
 
 
444
      if(.not.MA_Push_Get(MT_INT,ncont,'prim/shell',lprim,iprim))
 
445
     $     call errquit('wnbofile: lprim memory alloc. failed',0,
 
446
     &       MA_ERR)
 
447
 
 
448
      if(.not.MA_Push_Get(MT_INT,ncont+1,'contraction pointer',lptr,
 
449
     $     iptr))
 
450
     $     call errquit('wnbofile: lptr memory alloc. failed',0,
 
451
     &       MA_ERR)
 
452
 
 
453
      if(.not.MA_Push_Get(MT_INT,ncont,'contraction type',ltype,
 
454
     $     itype))
 
455
     $     call errquit('wnbofile: ltype memory alloc. failed',0,
 
456
     &       MA_ERR)
 
457
 
 
458
      if(.not.MA_Push_Get(MT_INT,ncont,'contraction/shell size',lncomp,
 
459
     $     incomp))
 
460
     $     call errquit('wnbofile: lncomp memory alloc. failed',0,
 
461
     &       MA_ERR)
 
462
 
 
463
      if(ga_nodeid().eq.0)then
 
464
        write(lfn,'(a)')'$BASIS'
 
465
        if(.not.MA_Alloc_Get(MT_Int,nbf,'scratch',lscr,iscr))
 
466
     $       call errquit('wnbofile: lscr memory alloc. failed',0,
 
467
     &       MA_ERR)
 
468
        do i = 1, nmo
 
469
          if(.not.bas_bf2ce(lbasis,i,int_mb(iscr+i-1)))
 
470
     $         call errquit('wnbofile: bas_bf2ce failed',0, BASIS_ERR)
 
471
        enddo 
 
472
 
 
473
        write(lfn,'(1x,a)')'CENTER = '
 
474
        write(lfn,'(5x,10i6)') (int_mb(iscr+i-1),i=1,nbf)
 
475
        
 
476
        i = 0
 
477
        nprim = 0
 
478
 
 
479
 
 
480
        int_mb(iptr) = 1
 
481
        maxL = 0
 
482
        do icont = 1,ncont
 
483
          if(.not.bas_continfo(lbasis,icont,type,int_mb(iprim+icont-1),
 
484
     $         ngeno,sphcart))
 
485
     $         call errquit('wnbofile: bas_continfo failed',0,
 
486
     &       BASIS_ERR)
 
487
          INT_mb(itype+icont-1) = type
 
488
          maxL = max(maxL,type)
 
489
          nprim = nprim + int_mb(iprim+icont-1)
 
490
          int_mb(iptr+icont) = 
 
491
     $         int_mb(iptr+icont-1) + int_mb(iprim+icont-1)
 
492
          index = iscr+i
 
493
          if(type.eq.-2)then
 
494
c
 
495
c>>>        spd shell
 
496
c
 
497
            int_mb(index)   = 1
 
498
            int_mb(index+1) = 101
 
499
            int_mb(index+2) = 102
 
500
            int_mb(index+3) = 103
 
501
 
 
502
            if( sphcart.eq.0)then
 
503
              int_mb(index+4) = 201
 
504
              int_mb(index+5) = 202
 
505
              int_mb(index+6) = 203
 
506
              int_mb(index+7) = 204
 
507
              int_mb(index+8) = 205
 
508
              int_mb(index+9) = 206
 
509
              i = i+10
 
510
              int_mb(incomp+icont-1) = 10
 
511
            elseif(sphcart.eq.1)then
 
512
              int_mb(index+4) = 251
 
513
              int_mb(index+5) = 252
 
514
              int_mb(index+6) = 253
 
515
              int_mb(index+7) = 254
 
516
              int_mb(index+8) = 255
 
517
              i = i+9
 
518
              int_mb(incomp+icont-1) = 10
 
519
            endif 
 
520
 
 
521
          elseif(type.eq.-1) then
 
522
c
 
523
c>>>        sp shell
 
524
c
 
525
            int_mb(index)   = 1
 
526
            int_mb(index+1) = 101
 
527
            int_mb(index+2) = 102
 
528
            int_mb(index+3) = 103
 
529
            i = i+4
 
530
            int_mb(incomp+icont-1) = 4
 
531
 
 
532
          elseif(type.eq.0) then
 
533
c
 
534
c>>>        s shell
 
535
c
 
536
            int_mb(index)   = 1
 
537
            i = i+1
 
538
            int_mb(incomp+icont-1) = 1
 
539
 
 
540
          elseif(type.eq.1) then
 
541
c
 
542
c>>>        p shell
 
543
c
 
544
            do j = 0,2
 
545
              int_mb(index+j) = 101+j
 
546
            enddo 
 
547
            i = i+3
 
548
            int_mb(incomp+icont-1) = 3
 
549
 
 
550
          elseif(type.eq.2) then
 
551
c
 
552
c>>>        d shell
 
553
c
 
554
            if( sphcart.eq.0)then
 
555
              do j = 0,5
 
556
                int_mb(index+j) = 201+j
 
557
              enddo 
 
558
              i = i+6
 
559
              int_mb(incomp+icont-1) = 6
 
560
            elseif(sphcart.eq.1)then
 
561
 
 
562
              do j = 0,4
 
563
                int_mb(index+j) = 251+j
 
564
              enddo 
 
565
              i = i+5
 
566
              int_mb(incomp+icont-1) = 5
 
567
            endif 
 
568
 
 
569
          elseif(type.eq.3) then
 
570
c
 
571
c>>>        f shell
 
572
c
 
573
            if(sphcart.eq.0)then
 
574
              do j = 0,9
 
575
                int_mb(index+j) = 301+j
 
576
              enddo 
 
577
              i = i+10
 
578
              int_mb(incomp+icont-1) = 10
 
579
            elseif(sphcart.eq.1)then
 
580
 
 
581
              do j = 0,6
 
582
                int_mb(index+j) = 351+j
 
583
              enddo 
 
584
              i = i+7
 
585
              int_mb(incomp+icont-1) = 7
 
586
            endif 
 
587
 
 
588
          elseif(type.eq.4) then
 
589
c
 
590
c>>>        g shell
 
591
c
 
592
            if(sphcart.eq.0)then
 
593
              do j = 0,14
 
594
                int_mb(index+j) = 401+j
 
595
              enddo 
 
596
              i = i+15
 
597
              int_mb(incomp+icont-1) = 15
 
598
            elseif(sphcart.eq.1)then
 
599
              do j = 0,8
 
600
                int_mb(index+j) = 451+j
 
601
              enddo
 
602
              i = i+9
 
603
              int_mb(incomp+icont-1) = 9
 
604
            endif 
 
605
 
 
606
          elseif(type.gt.4) then
 
607
            write(6,*)' only up to g functions allowed '
 
608
            call errquit('wnbofile: max angular momentum exceeded',0,
 
609
     &       BASIS_ERR)
 
610
          endif 
 
611
          
 
612
        enddo 
 
613
c
 
614
        write(lfn,'(2x,a)')'LABEL = '
 
615
        write(lfn,'(5x,10i6)') (int_mb(iscr+i-1),i=1,nbf)
 
616
        write(lfn,9200)
 
617
 
 
618
        if (.not. ma_free_heap(lscr))
 
619
     $       call errquit('wnbofile: lscr free heap failed',0, MA_ERR)
 
620
      endif
 
621
 
 
622
c
 
623
c>>>  End $BASIS keylist.
 
624
c
 
625
c
 
626
c>>> Begin $CONTRACT keylist.
 
627
      if (ga_nodeid().eq.0) then
 
628
        if(.not.MA_Alloc_Get(MT_dbl,nprim,'exponents',lexp,iexp))
 
629
     $       call errquit('wnbofile: lexp memory alloc. failed',0,
 
630
     &       MA_ERR)
 
631
        if(.not.MA_Alloc_Get(MT_dbl,nprim,'coefs',lcoef,icoef))
 
632
     $       call errquit('wnbofile: lcoef memory alloc failed',0,
 
633
     &       MA_ERR)
 
634
 
 
635
        index = 0
 
636
        do icont = 1,ncont
 
637
          if(.not.bas_get_exponent(lbasis,icont,dbl_mb(iexp+index)))
 
638
     $         call errquit('wnbofile: failed bas_get_exp',0, BASIS_ERR)
 
639
          if(.not.bas_get_coeff(lbasis,icont,dbl_mb(icoef+index)))
 
640
     $         call errquit('wnbofile: failed bas_get_coeff',0,
 
641
     &       BASIS_ERR)
 
642
          index = index + int_mb(iprim+icont-1)
 
643
        enddo 
 
644
 
 
645
        write(lfn,'(a)') '$CONTRACT'
 
646
        write(lfn,'(1x,a10,i5)') 'NSHELLS = ', ncont
 
647
        write(lfn,'(1x,a10,i5)') '   NEXP = ', nprim
 
648
        write(lfn,'(1x,a10)') '  NCOMP = '
 
649
        write(lfn,'(5x,10i6)') (int_mb(incomp+i),i=0,(ncont-1))
 
650
        write(lfn,'(1x,a10)') '  NPRIM = '
 
651
        write(lfn,'(5x,10i6)') (int_mb(iprim+i-1),i=1,ncont)
 
652
        write(lfn,'(1x,a10)') '   NPTR = '
 
653
        write(lfn,'(5x,10i6)') (int_mb(iptr+i-1),i=1,ncont)
 
654
        write(lfn,'(1x,a10)') '    EXP = '
 
655
        write(lfn,'(5x,4f20.10)') (dbl_mb(iexp+i-1),i=1,nprim)
 
656
 
 
657
        do l = 0,maxL
 
658
          write(lfn,'(5x,a)') cname(l+1)
 
659
c
 
660
c>>>      Reuse array space pointed to by iexp to print out
 
661
c>>>      s, p, d, f and g coefficients.
 
662
c
 
663
          call dcopy(nprim,0.0d0,0,dbl_mb(iexp),1)
 
664
          index = 0
 
665
          do icont = 1,ncont
 
666
            
 
667
            type=int_mb(itype+icont-1)
 
668
 
 
669
            sshell = (l.eq.0).and.(type.le.0)
 
670
            pshell = (l.eq.1).and.((type.eq.1).or.(type.le.-1))
 
671
            dshell = (l.eq.2).and.((type.eq.2).or.(type.le.-2))
 
672
            fshell = (l.eq.3).and.(type.eq.3)
 
673
            gshell = (l.eq.4).and.(type.eq.4)
 
674
 
 
675
            ilo = int_mb(iptr+icont-1)
 
676
            ihi = int_mb(iptr+icont)-1
 
677
c
 
678
c>>>        S-type coefficients.
 
679
c
 
680
            if(sshell)then
 
681
              do i = ilo,ihi
 
682
                dbl_mb(iexp+i-1)=dbl_mb(icoef+i-1)
 
683
              enddo 
 
684
            endif 
 
685
c
 
686
c>>>        P-type coefficients.
 
687
c
 
688
            if(pshell)then
 
689
              do i = ilo,ihi
 
690
                dbl_mb(iexp+i-1)=dbl_mb(icoef+i-1)
 
691
              enddo 
 
692
            endif 
 
693
c
 
694
c>>>        D-type coefficients.
 
695
c
 
696
            if(dshell)then
 
697
              do i = ilo,ihi
 
698
                dbl_mb(iexp+i-1)=dbl_mb(icoef+i-1)
 
699
              enddo 
 
700
            endif 
 
701
c
 
702
c>>>        F-type coefficients.
 
703
c
 
704
            if(fshell)then
 
705
              do i = ilo,ihi
 
706
                dbl_mb(iexp+i-1)=dbl_mb(icoef+i-1)
 
707
              enddo 
 
708
            endif 
 
709
c
 
710
c>>>        G-type coefficients.
 
711
c
 
712
            if(gshell)then
 
713
              do i = ilo,ihi
 
714
                dbl_mb(iexp+i-1)=dbl_mb(icoef+i-1)
 
715
              enddo 
 
716
            endif 
 
717
 
 
718
          enddo 
 
719
          
 
720
          write(lfn,'(5x,4f20.10)') (dbl_mb(iexp+i-1),i=1,nprim)
 
721
        enddo 
 
722
 
 
723
        write(lfn,9200)
 
724
 
 
725
        if (.not. ma_free_heap(lcoef))
 
726
     $       call errquit('wnbofile: lcoef free heap failed',0, MA_ERR)
 
727
 
 
728
        if (.not. ma_free_heap(lexp))
 
729
     $       call errquit('wnbofile: lexp free heap failed',0, MA_ERR)
 
730
      endif 
 
731
c
 
732
c>>> End $CONTRACT keylist.
 
733
c
 
734
c>>>  Read in molecular orbital vectors and occupation numbers.
 
735
c
 
736
      lenocc = nset*nmo
 
737
      if (.not. ma_push_get(mt_dbl, lenocc, 'wnbofile: mo evals',
 
738
     $     leval, ieval)) call errquit
 
739
     $     ('wnbofile: insufficient memory?', lenocc, MA_ERR)
 
740
 
 
741
      if (.not. ma_push_get(mt_dbl, lenocc, 'wnbofile: mo occ',
 
742
     $     locc, iocc)) call errquit
 
743
     $     ('wnbofile: mo occ insufficient memory?', lenocc, MA_ERR)
 
744
 
 
745
      if (.not. movecs_read(movecs_in, 1, dbl_mb(iocc), dbl_mb(ieval),
 
746
     $     g_movecs))
 
747
     $     call errquit('scf_movecs_read failed',0, DISK_ERR)
 
748
 
 
749
      if(.not.restricted)then
 
750
        if (.not. movecs_read(movecs_in, 2,
 
751
     $       dbl_mb(iocc+nbf), dbl_mb(ieval+nbf),
 
752
     $       g_movecs(2))) then
 
753
          call ga_copy(g_movecs(1), g_movecs(2))
 
754
          call dcopy(nbf,dbl_mb(iocc),1,dbl_mb(iocc+nbf),1)
 
755
          call dcopy(nbf,dbl_mb(ieval),1,dbl_mb(ieval+nbf),1)
 
756
        endif
 
757
      endif
 
758
c
 
759
c>>>    Begin $LCAOMO keylist.
 
760
c       
 
761
      if( ga_nodeid().eq.0)then
 
762
        if(.not.MA_Alloc_Get(MT_Dbl,nbf,'scratch',lscr,iscr))
 
763
     $       call errquit('wnbofile: lscr memory alloc. failed',0,
 
764
     &       MA_ERR)
 
765
 
 
766
        write(lfn,'(a)')'$LCAOMO'
 
767
 
 
768
        do icol = 1,nbf
 
769
          call ga_get(g_movecs(1),1,nmo,icol,icol,Dbl_MB(iscr),1)
 
770
          write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,nbf-1)
 
771
        enddo 
 
772
 
 
773
 
 
774
        if(.not. restricted) then
 
775
          do icol = 1,nbf
 
776
            call ga_get(g_movecs(2),1,nmo,icol,icol,Dbl_MB(iscr),1)
 
777
            write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,nbf-1)
 
778
          enddo
 
779
        endif
 
780
 
 
781
        write(lfn,9200)
 
782
        if (.not. ma_free_heap(lscr))
 
783
     $       call errquit('wnbofile: lscr free heap failed',0, MA_ERR)
 
784
      endif
 
785
c
 
786
c>>>  End LCAOMO keylist.
 
787
c
 
788
c
 
789
c>>>  Create bond order matrix from movecs array for density keylist.
 
790
c
 
791
c                     nmo   noc       t
 
792
c       rho(r,r') = N*SUM   SUM   C  C   chi(r) chi(r')
 
793
c                      jk    i     ji ik
 
794
c
 
795
c                                      /  noc     t  \
 
796
c       bond order matrix, gamma   =  | N*SUM C  C    |
 
797
c                               ij     \   i   ji ik /
 
798
c
 
799
      if( restricted) then
 
800
*ga:1:0
 
801
        if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: bond order', 
 
802
     $     0, 0, g_bo(1))) call errquit('wnbofile: g_bo(1)', 0, GA_ERR)
 
803
        call ga_zero(g_bo(1))
 
804
      else 
 
805
 
 
806
*ga:1:0
 
807
        if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: a bond order', 
 
808
     $     0, 0, g_bo(1))) call errquit('wnbofile: g_bo(1)', 0, GA_ERR)
 
809
        call ga_zero(g_bo(1))
 
810
        if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: b bond order', 
 
811
     $     0, 0, g_bo(2))) call errquit('wnbofile: g_bo(2)', 0, GA_ERR)
 
812
        call ga_zero(g_bo(2))
 
813
 
 
814
      endif 
 
815
c
 
816
c>>>  Create the density matrix in the MO representation, 
 
817
c>>>  a square matrix with orbital occupation 
 
818
c>>>  numbers on the diagonals and 0 elsewhere.  Store the matrix
 
819
c>>>  temporarily in the matrix reserved for the bond order matrix.
 
820
c
 
821
 
 
822
      if(.not.MA_push_get(MT_INT, nmo, 'indices 1', lind1, iind1))
 
823
     $     call errquit('wnbofile: lind1',0, MA_ERR)
 
824
      if(.not.MA_push_get(MT_INT, nmo, 'indices 2', lind2, iind2))
 
825
     $     call errquit('wnbofile: lind2',0, MA_ERR)
 
826
 
 
827
      if (ga_nodeid().eq.0)then
 
828
 
 
829
        do i=1, nmo
 
830
          int_mb( iind1+i-1 ) = i
 
831
          int_mb( iind2+i-1 ) = i
 
832
        enddo
 
833
 
 
834
*        write(6,*)' g_bo B4 zero'
 
835
*        call ga_print(g_bo(1))
 
836
*        call ga_zero (g_bo(1))
 
837
*        write(6,*)' g_bo after zero/B4 scatter'
 
838
*        call ga_print(g_bo(1))
 
839
        call ga_scatter ( g_bo(1), dbl_mb(iocc), int_mb(iind1),
 
840
     $       int_mb(iind2), nmo )
 
841
 
 
842
*        write(6,*)' g_bo after scatter'
 
843
*        call ga_print(g_bo(1))
 
844
      endif
 
845
 
 
846
*ga:1:0
 
847
      if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: scratch', 
 
848
     $     0, 0, g_scr)) call errquit('wnbofile: scratch', 0, GA_ERR)
 
849
      call ga_zero(g_scr)
 
850
c                                                t
 
851
c>>>  Form the product (Movecs)(Gamma_mo)(Movecs)  = Bond order matrix.
 
852
c
 
853
      call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_movecs(1), 
 
854
     $     g_bo(1), 0.0d0, g_scr)
 
855
      call ga_dgemm('n', 't', nmo, nmo, nmo, 1.0d0, g_scr, 
 
856
     $     g_movecs(1), 0.0d0, g_bo(1))
 
857
 
 
858
      if(.not. restricted) then
 
859
 
 
860
        if (ga_nodeid().eq.0)then
 
861
 
 
862
          do i=1, nmo
 
863
            int_mb( iind1+i-1 ) = i
 
864
            int_mb( iind2+i-1 ) = i
 
865
          enddo
 
866
 
 
867
          call ga_scatter ( g_bo(2), dbl_mb(iocc+nmo), int_mb(iind1),
 
868
     $         int_mb(iind2), nmo )
 
869
 
 
870
        endif
 
871
 
 
872
        call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_movecs(2), 
 
873
     $       g_bo(2), 0.0d0, g_scr)
 
874
        call ga_dgemm('n', 't', nmo, nmo, nmo, 1.0d0, g_scr, 
 
875
     $       g_movecs(2), 0.0d0, g_bo(2))
 
876
      endif 
 
877
c
 
878
c>>>  Begin $DENSITY keylist.
 
879
c
 
880
      if( ga_nodeid().eq.0)then
 
881
        if(.not.MA_Alloc_Get(MT_Dbl,nbf,'scratch',lscr,iscr))
 
882
     $       call errquit('wnbofile: lscr memory alloc. failed',0,
 
883
     &       MA_ERR)
 
884
 
 
885
        write(lfn,'(a)')'$DENSITY'
 
886
c
 
887
c>>>  Write out lower triangular section of alpha density matrix.
 
888
c
 
889
 
 
890
        do icol = 1,nbf
 
891
          call ga_get(g_bo(1),1,nmo,icol,icol,Dbl_MB(iscr),1)
 
892
          write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,icol-1)
 
893
        enddo 
 
894
 
 
895
 
 
896
        if(.not. restricted) then
 
897
          do icol = 1,nbf
 
898
            call ga_get(g_bo(2),1,nmo,icol,icol,Dbl_MB(iscr),1)
 
899
            write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,icol-1)
 
900
          enddo
 
901
        endif
 
902
 
 
903
        write(lfn,9200)
 
904
        if (.not. ma_free_heap(lscr))
 
905
     $       call errquit('wnbofile: lscr free heap failed',0, MA_ERR)
 
906
      endif
 
907
c
 
908
c>>>  End $DENSITY keylist.
 
909
c
 
910
c
 
911
c>>>  Create overlap matrix.
 
912
c
 
913
      call int_init(rtdb,1,lbasis)
 
914
      if (.not.int_normalize(rtdb,lbasis)) call errquit
 
915
     &    ('wnbofile: int_normalize failed',911, INT_ERR)
 
916
      g_over  = ga_create_atom_blocked(lgeom, lbasis,'Overlap')
 
917
c
 
918
      call ga_zero(g_over)
 
919
      call int_1e_ga(lbasis, lbasis, g_over, 'overlap', .false.)
 
920
c
 
921
c>>>  Begin $OVERLAP keylist.
 
922
c
 
923
      if( ga_nodeid().eq.0)then
 
924
        if(.not.MA_Alloc_Get(MT_Dbl,nbf,'scratch',lscr,iscr))
 
925
     $       call errquit('wnbofile: lscr memory alloc. failed',0,
 
926
     &       MA_ERR)
 
927
 
 
928
        write(lfn,'(a)')'$OVERLAP'
 
929
c
 
930
c>>>  Write out lower triangular section of overlap matrix.
 
931
c
 
932
        do icol = 1,nbf
 
933
          call ga_get(g_over,1,nmo,icol,icol,Dbl_MB(iscr),1)
 
934
          write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,icol-1)
 
935
        enddo 
 
936
 
 
937
        write(lfn,9200)
 
938
        if (.not. ma_free_heap(lscr))
 
939
     $       call errquit('wnbofile: lscr free heap failed',0, MA_ERR)
 
940
 
 
941
      endif 
 
942
c
 
943
c>>>  End $OVERLAP keylist.
 
944
c
 
945
c
 
946
c>>>  Create Fock matrix in AO basis.
 
947
c
 
948
      if(restricted) then
 
949
*ga:1:0
 
950
        if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: Fock matrix', 
 
951
     $     0, 0, g_fock(1))) call errquit('wnbofile: g_fock(1)', 0,
 
952
     &       GA_ERR)
 
953
        call ga_zero(g_fock(1))
 
954
      else 
 
955
 
 
956
*ga:1:0
 
957
        if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: alpha Fock', 
 
958
     $     0, 0, g_fock(1))) call errquit('wnbofile: g_fock(1)', 0,
 
959
     &       GA_ERR)
 
960
        call ga_zero(g_fock(1))
 
961
*ga:1:0
 
962
        if (.not. ga_create(MT_DBL, nmo, nmo, 'wnbofile: beta Fock', 
 
963
     $     0, 0, g_fock(2))) call errquit('wnbofile: g_fock(2)', 0,
 
964
     &       GA_ERR)
 
965
        call ga_zero(g_fock(2))
 
966
 
 
967
      endif 
 
968
c
 
969
c>>>  Create Fock matrix in MO representation, 
 
970
c>>>  a square matrix with oribital energies
 
971
c>>>  on the diagonals and 0 elsewhere.
 
972
c
 
973
      if (ga_nodeid().eq.0)then
 
974
 
 
975
        do i=1, nmo
 
976
          int_mb( iind1+i-1 ) = i
 
977
          int_mb( iind2+i-1 ) = i
 
978
        enddo
 
979
 
 
980
        call ga_scatter ( g_fock(1), dbl_mb(ieval), int_mb(iind1),
 
981
     $       int_mb(iind2), nmo )
 
982
 
 
983
      endif
 
984
c
 
985
c>>>  Create Fock Matrix in AO representation.
 
986
c                                   t
 
987
c>>>  Fock_ao = S*Movecs*Fock_mo*Movecs *S
 
988
c
 
989
      call ga_dgemm('n', 't', nmo, nmo, nmo, 1.0d0, g_fock(1), 
 
990
     $       g_movecs(1), 0.0d0, g_scr)
 
991
      call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_movecs(1), 
 
992
     $       g_scr, 0.0d0, g_fock(1))
 
993
      call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_fock(1), 
 
994
     $       g_over, 0.0d0, g_scr)
 
995
      call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_over, 
 
996
     $       g_scr, 0.0d0, g_fock(1))
 
997
 
 
998
      if(.not. restricted) then
 
999
        
 
1000
        if (ga_nodeid().eq.0)then
 
1001
 
 
1002
          do i=1, nmo
 
1003
            int_mb( iind1+i-1 ) = i
 
1004
            int_mb( iind2+i-1 ) = i
 
1005
          enddo
 
1006
 
 
1007
          call ga_scatter ( g_fock(2), dbl_mb(ieval+nmo), int_mb(iind1),
 
1008
     $         int_mb(iind2), nmo )
 
1009
 
 
1010
        endif
 
1011
 
 
1012
        call ga_dgemm('n', 't', nmo, nmo, nmo, 1.0d0, g_fock(2), 
 
1013
     $       g_movecs(2), 0.0d0, g_scr)
 
1014
        call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_movecs(2), 
 
1015
     $       g_scr, 0.0d0, g_fock(2))
 
1016
        call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_fock(2), 
 
1017
     $       g_over, 0.0d0, g_scr)
 
1018
        call ga_dgemm('n', 'n', nmo, nmo, nmo, 1.0d0, g_over, 
 
1019
     $       g_scr, 0.0d0, g_fock(2))
 
1020
 
 
1021
      endif 
 
1022
c
 
1023
c>>>  Begin $FOCK keylist.
 
1024
c
 
1025
      if( ga_nodeid().eq.0)then
 
1026
        if(.not.MA_Alloc_Get(MT_Dbl,nbf,'scratch',lscr,iscr))
 
1027
     $       call errquit('wnbofile: lscr memory alloc. failed',0,
 
1028
     &       MA_ERR)
 
1029
 
 
1030
        write(lfn,'(a)')'$FOCK'
 
1031
c
 
1032
c>>>  Write out lower triangular section of AO Fock matrix.
 
1033
c
 
1034
        do icol = 1,nbf
 
1035
          call ga_get(g_fock(1),1,nmo,icol,icol,Dbl_MB(iscr),1)
 
1036
          write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,icol-1)
 
1037
        enddo 
 
1038
 
 
1039
        if(.not. restricted) then
 
1040
          do icol = 1,nbf
 
1041
            call ga_get(g_fock(2),1,nmo,icol,icol,Dbl_MB(iscr),1)
 
1042
            write(lfn,'(5f20.10)')(Dbl_MB(iscr+i),i=0,icol-1)
 
1043
          enddo
 
1044
        endif
 
1045
 
 
1046
        write(lfn,9200)
 
1047
        if (.not. ma_free_heap(lscr))
 
1048
     $       call errquit('wnbofile: lscr free heap failed',0, MA_ERR)
 
1049
      endif
 
1050
 
 
1051
c
 
1052
c>>>  End $FOCK keylist.
 
1053
c
 
1054
c
 
1055
c>>>  Print this out for debugging purposes.
 
1056
 
 
1057
      close(lfn)
 
1058
c
 
1059
c>>>  Free handles and destroy arrays.
 
1060
c
 
1061
      if (restricted) then
 
1062
        if (.not. ga_destroy(g_scr)) 
 
1063
     $       call errquit('wnbofile: destroy g_scr',0, GA_ERR)
 
1064
        if (.not. ga_destroy(g_movecs(1))) 
 
1065
     $       call errquit('wnbofile: destroy g_movecs',0, GA_ERR)
 
1066
        if (.not. ga_destroy(g_bo(1))) 
 
1067
     $       call errquit('wnbofile: destroy g_bo',0, GA_ERR)
 
1068
        if (.not. ga_destroy(g_over)) 
 
1069
     $     call errquit('wnbofile: destroy g_over',0, GA_ERR)
 
1070
        if (.not. ga_destroy(g_fock(1))) 
 
1071
     $       call errquit('wnbofile: destroy g_fock',0, GA_ERR)
 
1072
      else
 
1073
        if (.not. ga_destroy(g_scr)) 
 
1074
     $       call errquit('wnbofile: destroy g_scr',0, GA_ERR)
 
1075
        if (.not. ga_destroy(g_movecs(2)))
 
1076
     $       call errquit('wnbofile: destroy beta MOs', 0, GA_ERR)
 
1077
        if (.not. ga_destroy(g_movecs(1)))
 
1078
     $       call errquit('wnbofile: destroy alpha MOs', 0, GA_ERR)
 
1079
        if (.not. ga_destroy(g_over)) 
 
1080
     $     call errquit('wnbofile: destroy g_over',0, GA_ERR)
 
1081
        if (.not. ga_destroy(g_bo(2)))
 
1082
     $       call errquit('wnbofile: destroy beta BO matrix', 0, GA_ERR)
 
1083
        if (.not. ga_destroy(g_bo(1)))
 
1084
     $       call errquit('wnbofile: destroy alpha BO matrix', 0,
 
1085
     &       GA_ERR)
 
1086
        if (.not. ga_destroy(g_fock(1))) 
 
1087
     $       call errquit('wnbofile: destroy g_fock',0, GA_ERR)
 
1088
        if (.not. ga_destroy(g_fock(2))) 
 
1089
     $       call errquit('wnbofile: destroy g_fock',0, GA_ERR)
 
1090
      endif
 
1091
 
 
1092
      if (.not. ma_pop_stack(lind2))
 
1093
     $     call errquit('wnbofile: lind2 pop stack failed',0, MA_ERR)
 
1094
 
 
1095
      if (.not. ma_pop_stack(lind1))
 
1096
     $     call errquit('wnbofile: lind1 pop stack failed',0, MA_ERR)
 
1097
 
 
1098
      if (.not. ma_pop_stack(locc))
 
1099
     $     call errquit('wnbofile: locc pop stack failed',0, MA_ERR)
 
1100
 
 
1101
      if (.not. ma_pop_stack(leval))
 
1102
     $     call errquit('wnbofile: leval pop stack failed',0, MA_ERR)
 
1103
 
 
1104
      if (.not. ma_pop_stack(lncomp))
 
1105
     $     call errquit('wnbofile: lncomp pop stack failed',0, MA_ERR)
 
1106
 
 
1107
      if (.not. ma_pop_stack(ltype))
 
1108
     $     call errquit('wnbofile: ltype pop stack failed',0, MA_ERR)
 
1109
 
 
1110
      if (.not. ma_pop_stack(lptr))
 
1111
     $     call errquit('wnbofile: lptr pop stack failed',0, MA_ERR)
 
1112
 
 
1113
      if (.not. ma_pop_stack(lprim))
 
1114
     $     call errquit('wnbofile: lprim pop stack failed',0, MA_ERR)
 
1115
 
 
1116
      if (.not. ma_pop_stack(ltags))
 
1117
     $     call errquit('wnbofile: ltags pop stack failed',0, MA_ERR)
 
1118
 
 
1119
      if (.not. ma_pop_stack(lchg))
 
1120
     $     call errquit('wnbofile: lchg pop stack failed',0, MA_ERR)
 
1121
 
 
1122
      if (.not. ma_pop_stack(lcoord))
 
1123
     $     call errquit('wnbofile: lcoord pop stack failed',0, MA_ERR)
 
1124
c ------  Deallocate memory ----------
 
1125
      call int_terminate()
 
1126
   
 
1127
 147  continue
 
1128
 
 
1129
      if (.not. bas_destroy(lbasis)) call errquit
 
1130
     $     ('wnbofile: basis destroy failed',0, BASIS_ERR)
 
1131
 
 
1132
      if (.not. geom_destroy(lgeom)) call errquit
 
1133
     $     ('wnbofile: geom destroy failed', 0, GEOM_ERR)
 
1134
 
 
1135
      call ga_sync()
 
1136
 
 
1137
 9200 format('$END')
 
1138
 9250 format(/,1x,'Input for gennbo program written to file ',a,'.')
 
1139
      return       
 
1140
      end