1      subroutine rdexbas(basfile)
2      implicit double precision (a-h,o-z)
3      parameter (maxbas=15)
4      parameter (mxel=100)
5      parameter (mxaorb=40)
6      common /basext/ ibasin(mxel),eborbs(maxbas),
7     &     pbas(mxaorb,mxaorb,maxbas)
8      common /basopt/ extbas
9      logical extbas
10      character*75 basfile,line
11      character*8 aname
12
13      extbas = .false.
14
15      open (unit=99,file=basfile,status='old',err=90)
16      num = 1
17      do while(.true.)
18
19         read (99,'(A8)',end=10,err=100) aname
20         read (99,*,err=100) line
21         write(6,*) aname
22         read (99,'(i4,i10)',err=100) icharge,norbs
23         ibasin(icharge) = num
24         eborbs(num) = norbs
25         write(6,*) 'norbs=', eborbs(num),'icharge=', icharge
26         do i=1,norbs
27            read (99,'(6f15.10)',err=100) (pbas(i,j,num),j=1,norbs)
28         end do
29         num = num + 1
30         read(99,*,err=100) line
31      end do
32      close(99)
33
3410    continue
35      if (num.gt.1) extbas = .true.
36      return
37
3890    continue
39      print*,'error opening file'
40      return
41
42100   continue
43      return
44      end
45
46      subroutine newdenmad(
47     &                     p,q)
48c THIS IS REALLY newdenmat
49
50c     returns density matrix in p common /densty/
51      implicit double precision (a-h,o-z)
52      parameter (numatm=2000)
53      parameter (maxbas=15)
54      parameter (mxel=100)
55      parameter (mxaorb=40)
56      logical valenc,bonds,ovrlap,atomic,doori,orient,dolap,doelf
57      common /rdwr/   iun1,iun2,iun3,iun4,iun5
58      common /moldat/ natoms, norbs, nelecs,nat(numatm)
59      common /option/ ipsi,idebug,valenc,bonds,ovrlap,atomic,doori,
60     &                dolap,doelf
61
62      common /basext/ ibasin(mxel),eborbs(maxbas),
63     &     pbas(mxaorb,mxaorb,maxbas)
64      common /orbhlp/ mxorb,iuhf,ispd
65
66      character*2 elemnt
67      common /elem/elemnt(mxel)
68      dimension p(*),q(*)
69
70      do i=1,norbs
71         do j=1,norbs
72            q((j-1)*mxorb + i) = 0.0d0
73         end do
74      end do
75
76      ix = 0
77      do 700 k=1,natoms
78         n = eborbs(ibasin(nat(k)))
79         write(6,*) 'n=', n, 'K=',k
80         if (ovrlap) then
81c
82c============ do overlap densities only ========================
83c
84            do i=1,n
85               do j=1,n
86                     p((ix+j-1)*mxorb + (ix+i)) = 0.0d0
87                  end do
88               end do
89            else
90c
91c============ molecule minus atoms =============================
92c
93               orient = .false.
94C               do i=1,norien
95C                  if (ori(i).eq.k) then
96C                     orient = .true.
97C                     ibal   = i
98C                  endif
99C               end do
100               if (orient.and.doori) then
101                  write (6,*) 'Orienting '
102C
103               else
104c
105c----- for all other atoms substract spherical density -----
106c
107                  do i=1,n
108                     do j=1,n
109c                        write (6,*) 'before:', ix,i,j,p((ix+j-1)*mxorb
110c     $                  + (ix+i)), pbas((ix+i),(ix+j),ibasin(nat(k)))
111                        p((ix+j-1)*mxorb + (ix+i)) =
112     $                   p((ix+j-1)*mxorb + (ix+i)) - pbas((ix+i),(ix+j)
113     $                                 ,ibasin(nat(k)))
114                        q((ix+j-1)*mxorb + (ix+i)) =
115     $                        pbas((ix+i),(ix+j),ibasin(nat(k)))
116C                        write (6,*) 'after:' ,p(ix+i,ix+j)
117                 end do
118              end do
119           endif
120c
121c============== clear out the overlap region of the density matrix
122c               if ATOMIC is supplied
123c
124           if (atomic) then
125              idavje = n
126              do i=ix+1,ix+idavje
127                 if (k.ne.0) then
128                    do j=1,ix
129                       p((j-1)*mxorb + i) = 0.0d0
130                    end do
131                 endif
132                 if (k.ne.natoms) then
133                    do j=ix+idavje+1,norbs
134                       p((j-1)*mxorb + i) = 0.0d0
135                    end do
136                 endif
137              end do
138           endif
139        endif
140        ix = ix + n
141 700  continue
142      if (bonds.and.idebug.eq.1) then
143         write(iun3,*)' '
144         write(iun3,*)'***** Atomic Density Matrix *****'
145         write(iun3,*)' '
146         call prev(q,norbs,norbs,mxorb)
147         write(iun3,*)' '
148      endif
149
150c############### end bonds , overlap , atomic ##################
151      if (idebug.eq.1) then
152          write(iun3,'(''   DENSITY MATRIX USED BY MAP'')')
153          call prev(p,norbs,norbs,mxorb)
154      endif
155
156      return
157      end
158