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