1      Subroutine dft_getspm(geom,lmax,finest,iga_dens,basis)
2
3C$Id$
4      implicit none
5#include "errquit.fh"
6      integer lmax               ! max value of ang momentum
7      integer iga_dens           ! GA handle for DM
8      integer basis              ! basis handle
9      integer finest             ! finest level for CMM
10      integer  geom
11C
12C     local
13C
14      double precision centerl(3),q
15      double precision total_nuc_charge
16      integer natoms,iatom
17      integer max_at_bf,max_at_bf2
18      integer lcoord,icoord,lcharge,icharge,itags,ltags
19      integer ldipole,idipole,lPmat,iPmat
20      integer lxcrd,ixcrd,lycrd,iycrd,lzcrd,izcrd
21      integer i,ihi,istart,ilo,icent
22
23#include "bas.fh"
24#include "geom.fh"
25
26#include "mafdecls.fh"
27#include "stdio.fh"
28
29
30      if(.not. geom_ncent(geom, natoms))
31     &  call errquit('dft_scf: geom_ncent failed',0, GEOM_ERR)
32
33      if(.not.Ma_Alloc_Get(MT_Dbl,natoms*3,'coordinates',lcoord,icoord))
34     &  call errquit('dft_getspm: cannot allocate coordinates',0,
35     &       MA_ERR)
36      if(.not.Ma_Push_Get(MT_Dbl,natoms,'charges',lcharge,icharge))
37     &  call errquit('dft_getspm: cannot allocate charges',0,
38     &       MA_ERR)
39      if(.not.Ma_Push_Get(MT_Byte,natoms*16,'center tags',ltags,itags))
40     &  call errquit('dft_getspm: cannot allocate center tags',0,
41     &       MA_ERR)
42
43      if(.not. geom_cart_get(geom, natoms, Byte_MB(itags),
44     &  Dbl_MB(icoord), Dbl_MB(icharge)))
45     &  call errquit('dft_getspm: geom_cart_get failed',1, GEOM_ERR)
46
47c****
48c**** compute center of charge
49c****
50
51      centerl(1) = 0.d0
52      centerl(2) = 0.d0
53      centerl(3) = 0.d0
54
55      do icent = 0,natoms-1
56        q =dbl_mb(icharge+icent)
57        centerl(1) = centerl(1) + q*dbl_mb(icoord+3*icent  )
58        centerl(2) = centerl(2) + q*dbl_mb(icoord+3*icent+1)
59        centerl(3) = centerl(3) + q*dbl_mb(icoord+3*icent+2)
60      enddo
61      if(.not. geom_nuc_charge(geom,total_nuc_charge)) call errquit
62     &  ('dft_getspm: geom_nuc_charge failed',20, GEOM_ERR)
63
64      centerl(1) = centerl(1)/total_nuc_charge
65      centerl(2) = centerl(2)/total_nuc_charge
66      centerl(3) = centerl(3)/total_nuc_charge
67
68      if(.not.ma_pop_stack(ltags))
69     &  call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
70
71      if(lmax.gt.0) then
72        if(.not.Ma_Push_Get(MT_Dbl,3*natoms,'dipole',ldipole,
73     &    idipole))
74     &    call errquit('dft_getspm: cannot allocate dipole',0, MA_ERR)
75      endif
76
77      max_at_bf = 0
78      do iatom = 1, natoms
79        if(.not. bas_ce2bfr(basis, iatom, ilo, ihi))
80     $    call errquit('quadvxc0: bas_ce2bfr failed', iatom, BASIS_ERR)
81        max_at_bf = max(max_at_bf, ihi-ilo+1)
82      enddo
83      max_at_bf2=max_at_bf*max_at_bf
84
85      if(.not.Ma_Push_Get(MT_Dbl,max_at_bf2,'local DM',lPmat,iPmat))
86     &  call errquit('dft_getspm: cannot allocate local DM',0, MA_ERR)
87
88      call dft_genspm(lmax,iga_dens,basis,
89     &  natoms,centerl,dbl_MB(iPmat),max_at_bf,
90     &  DBL_MB(icharge),DBL_MB(idipole),dbl_mb(icoord))
91c
92c     move the coordinates form coord(3,nat) to x y z arrays
93c
94      if(.not.Ma_Push_Get(MT_Dbl,natoms,'x coord',lxcrd,ixcrd))
95     &  call errquit('dft_getspm: cannot allocate x coord',0, MA_ERR)
96      if(.not.Ma_Push_Get(MT_Dbl,natoms,'y coord',lycrd,iycrd))
97     &  call errquit('dft_getspm: cannot allocate y coord',0, MA_ERR)
98      if(.not.Ma_Push_Get(MT_Dbl,natoms,'z coord',lzcrd,izcrd))
99     &  call errquit('dft_getspm: cannot allocate z coord',0, MA_ERR)
100      istart=0
101      do i=0,natoms-1
102        dbl_mb(ixcrd+i)=dbl_mb(icoord+istart)
103        dbl_mb(iycrd+i)=dbl_mb(icoord+istart+1)
104        dbl_mb(izcrd+i)=dbl_mb(icoord+istart+2)
105        istart=istart+3
106      enddo
107
108      if(.not.ma_free_heap(lcoord))
109     &  call errquit('dft_getspm: cannot free heap',0, MA_ERR)
110c
111c     re-center using geometrical center
112c
113      call dft_geocent(natoms,
114     &  dbl_mb(ixcrd),dbl_mb(iycrd),dbl_mb(izcrd))
115c
116c     call to CMM code
117c
118C      call cmm_cmm(natoms,dbl_mb(ixcrd),dbl_mb(iycrd),dbl_mb(izcrd),
119C     &  dbl_mb(icharge),finest)
120
121      if(.not.ma_pop_stack(lzcrd))
122     &  call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
123      if(.not.ma_pop_stack(lycrd))
124     &  call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
125      if(.not.ma_pop_stack(lxcrd))
126     &  call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
127      if(.not.ma_pop_stack(lPmat))
128     &  call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
129      if(lmax.gt.0) then
130        if(.not.ma_pop_stack(ldipole))
131     &    call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
132      endif
133      if(.not.ma_pop_stack(lcharge))
134     &    call errquit('dft_getspm: cannot pop stack',0, MA_ERR)
135c
136c       compute geometrical center for CMM
137c
138
139      return
140      end
141      subroutine dft_geocent(natoms,xcrd,ycrd,zcrd)
142      implicit none
143c-----------------------------------------------------------------------
144c     Find the ranges of the molecular system in the x, y, and z directions
145c-----------------------------------------------------------------------
146#include "stdio.fh"
147#include "global.fh"
148#include "tcgmsg.fh"
149      integer natoms
150      double precision xcrd(natoms),ycrd(natoms),zcrd(natoms)
151c
152      double precision xmin,xmax
153      double precision ymin,ymax
154      double precision zmin,zmax
155      double precision xorigin,yorigin,zorigin
156      double precision xrange,yrange,zrange
157      double precision two
158      parameter (two=2.d0)
159
160      integer i
161
162      xmin = 1.e+30
163      ymin = 1.e+30
164      zmin = 1.e+30
165      xmax = -1.e+30
166      ymax = -1.e+30
167      zmax = -1.e+30
168
169      do i = 1 , natoms
170
171        xmin = min(xmin,xcrd(i))
172        ymin = min(ymin,ycrd(i))
173        zmin = min(zmin,zcrd(i))
174        xmax = max(xmax,xcrd(i))
175        ymax = max(ymax,ycrd(i))
176        zmax = max(zmax,zcrd(i))
177
178      end do
179
180      xrange = dabs(xmin-xmax)
181      yrange = dabs(ymin-ymax)
182      zrange = dabs(zmin-zmax)
183
184c-----------------------------------------------------------------------
185c     Compute the geometrical center of molecular system
186c-----------------------------------------------------------------------
187
188      xorigin = (xmin + xmax)/TWO
189      yorigin = (ymin + ymax)/TWO
190      zorigin = (zmin + zmax)/TWO
191
192c-----------------------------------------------------------------------
193c     Refer the coordinates to the geometrical center of molecular system
194c-----------------------------------------------------------------------
195
196      do i = 1 , natoms
197
198        xcrd(i) =  xcrd(i) - xorigin
199        ycrd(i) =  ycrd(i) - yorigin
200        zcrd(i) =  zcrd(i) - zorigin
201
202      end do
203
204c-----------------------------------------------------------------------
205c     Output for checking
206c-----------------------------------------------------------------------
207      if(ga_nodeid().eq.0) then
208        write(LuOut, '("minima",3(1x,f8.3))') xmin, ymin, zmin
209        write(LuOut, '("maxima",3(1x,f8.3)/)') xmax, ymax, zmax
210
211        write(LuOut, '("xwidth",f8.3)')  xrange
212        write(LuOut, '("ywidth",f8.3)')  yrange
213        write(LuOut, '("zwidth",f8.3/)') zrange
214
215        write(LuOut, '("origin",3(f8.3,1x)//)')
216     &       xorigin, yorigin, zorigin
217      endif
218
219      return
220      end
221