1 subroutine hnd_mtpole(rtdb,basis,geom,mpole) 2c 3c $Id$ 4c 5c This routine calculates the multipole moments around 6c choosen center of expansion 7c 8c mpole = 1 : Dipole 9c mpole = 2 : Quadrupole 10c mpole = 3 : Octupole 11c 12 implicit none 13#include "errquit.fh" 14#include "global.fh" 15#include "mafdecls.fh" 16#include "nwc_const.fh" 17#include "stdio.fh" 18#include "geom.fh" 19#include "rtdb.fh" 20#include "cosmo.fh" 21c 22 integer rtdb ! [input] rtdb handle 23 integer basis ! [input] basis handle 24 integer geom ! [input] geometry handle 25 integer mpole ! [input] order of multipole 26c 27c Setting a max_pole. We have up to Octupole. 28c 29 integer max_pole, maxcomp 30 parameter (max_pole = 3, maxcomp = (max_pole+1)*(max_pole+2)/2) 31c 32 double precision cm(3), mptval(maxcomp,3) 33 double precision mptval2(maxcomp,3) 34 double precision octo, debye, buck, cx, cy, cz 35 double precision dipol, dipefc, diptot, rsquar, dum1, dum2 36 logical status 37 integer Nxyz(3) 38 integer ncomp, nat, iat, icomp, i,j 39 integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt 40 integer nefc, l_efcc, k_efcc, l_efcz, k_efcz 41 integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2) 42 character*3 scftyp 43 character*16 at_tag 44 data octo /0.711688d+00/ 45 data buck /1.344911d+00/ 46 data debye /2.54176568d+00/ 47c 48 if (ga_nodeid().eq.0.and.mpole.eq.1) write(luout,2001) 49 if (ga_nodeid().eq.0.and.mpole.eq.2) write(luout,3001) 50 if (ga_nodeid().eq.0.and.mpole.eq.3) write(luout,9001) 51c 52c Initialize integrals 53c 54 call int_init(rtdb,1, basis) 55 call schwarz_init(geom, basis) 56c 57c Get density matrix 58c 59 call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp, 60 & nclosed,nopen,nvirt) 61c 62c Determine number of components and zero mptval 63c 64 ncomp = (mpole+1)*(mpole+2)/2 65 call dcopy(2*maxcomp,0.0d0,0,mptval,1) 66c 67c Get center of expansion 68c 69 call prp_cent(rtdb,geom,cm) 70c 71c Electronic contribution 72c 73 call hnd_mtpcon(basis,geom,g_dens(ndens),mptval,mpole,cm) 74c 75c Only node 0 writing to rtdb 76c 77 status = rtdb_parallel(.false.) 78c 79 if (ga_nodeid().gt.0) goto 1000 80c 81c Allocate some local memory 82c 83 if (.not.geom_ncent(geom,nat)) call 84 & errquit('hnd_efgmap: geom_ncent failed',911,GEOM_ERR) 85 if (.not. ma_push_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 86 & call errquit('hnd_mtpole: ma failed',911,MA_ERR) 87 if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 88 & call errquit('hnd_mtpole: ma failed',911,MA_ERR) 89c 90 do iat=1,nat 91 if (.not.geom_cent_get(geom,iat,at_tag, 92 & dbl_mb(k_xyzpt+3*(iat-1)),dbl_mb(k_zanpt+iat-1))) call 93 & errquit('hnd_mtpole: geom_cent_get failed',911,GEOM_ERR) 94 enddo ! iat 95c 96c For quadrupole case, need electronic contribution for 97c diamagnetic susceptibility 98c 99 if (mpole.eq.2) rsquar=mptval(1,1)+mptval(4,1)+mptval(6,1) 100c 101c Nuclear contribution 102c 103 do icomp = 1, ncomp 104 mptval(icomp,1) = -mptval(icomp,1) 105 call getNxyz(mpole,icomp,Nxyz) 106 do iat = 1, nat 107 cx = 1d0 108 cy = 1d0 109 cz = 1d0 110 if(Nxyz(1).ne.0) 111 & cx = (dbl_mb(k_xyzpt +3*(iat-1))-cm(1))**Nxyz(1) 112 if(Nxyz(2).ne.0) 113 & cy = (dbl_mb(k_xyzpt+1+3*(iat-1))-cm(2))**Nxyz(2) 114 if(Nxyz(3).ne.0) 115 & cz = (dbl_mb(k_xyzpt+2+3*(iat-1))-cm(3))**Nxyz(3) 116 mptval(icomp,1) = mptval(icomp,1) + 117 & dbl_mb(k_zanpt+iat-1) * cx * cy * cz 118 enddo ! iat 119 enddo ! icomp 120c 121c ----- form -efc- contribution ----- 122c from cosmo point charges !!!! 123c 124 if (cosmo_last) then 125 if (.not.rtdb_get(rtdb,'cosmo:nefc',mt_int,1,nefc)) 126 & call errquit('hnd_mtpole: rtdb get failed for nefc ',911, 127 & rtdb_err) 128 if (.not.ma_push_get(mt_dbl,nefc*3,'efcc',l_efcc,k_efcc)) 129 & call errquit('hnd_mtpole: malloc k_efcc fail',911,ma_err) 130 if (.not.ma_push_get(mt_dbl,nefc,'efcz',l_efcz,k_efcz)) 131 & call errquit('hnd_mtpole: malloc k_efcz fail',911,ma_err) 132 if (.not.rtdb_get(rtdb,'cosmo:efcc',mt_dbl,3*nefc, 133 & dbl_mb(k_efcc))) call 134 & errquit('hnd_mtpole: rtdb get failed efcc',912,rtdb_err) 135 if (.not.rtdb_get(rtdb,'cosmo:efcz',mt_dbl,nefc, 136 & dbl_mb(k_efcz))) call 137 & errquit('hnd_mtpole: rtdb get failed efcz',913,rtdb_err) 138 do icomp = 1, ncomp 139 call getNxyz(2,icomp,Nxyz) 140 do iat=1,nefc 141 cx=1d0 142 cy=1d0 143 cz=1d0 144 if(Nxyz(1).ne.0) 145 & cx = (dbl_mb(k_efcc+3*(icomp-1) )-cm(1))**Nxyz(1) 146 if(Nxyz(2).ne.0) 147 & cy = (dbl_mb(k_efcc+3*(icomp-1)+1)-cm(2))**Nxyz(2) 148 if(Nxyz(3).ne.0) 149 & cz = (dbl_mb(k_efcc+3*(icomp-1)+2)-cm(3))**Nxyz(3) 150 mptval(icomp,2) = mptval(icomp,2) 151 & + dbl_mb(k_efcz+iat-1)*cx*cy*cz 152 enddo ! iat 153 enddo ! icomp 154 if (.not.ma_chop_stack(l_efcc)) call 155 & errquit('hnd_mtpole: chop stack l_efcc',913,ma_err) 156 endif 157c 158c Print output for 159c mpole = 1 : Dipole 160c mpole = 2 : Quadrupole 161c mpole = 3 : Octupole 162c 163 goto (101,102,103) mpole 164c 165c Dipole output 166c 167 101 dipol = sqrt(mptval(1,1)*mptval(1,1) + 168 & mptval(2,1)*mptval(2,1) + 169 & mptval(3,1)*mptval(3,1)) 170 dipefc = sqrt(mptval(1,2)*mptval(1,2) + 171 & mptval(2,2)*mptval(2,2) + 172 & mptval(3,2)*mptval(3,2)) 173 diptot = sqrt((mptval(1,1)+mptval(1,2))**2 + 174 & (mptval(2,1)+mptval(2,2))**2 + 175 & (mptval(3,1)+mptval(3,2))**2) 176 write(luout,2995) dipol,mptval(1,1),mptval(1,2),mptval(2,1), 177 & mptval(2,2),mptval(3,1),mptval(3,2),dipefc, 178 & diptot 179c 180 do icomp = 1, ncomp 181 mptval(icomp,1) = mptval(icomp,1) * debye 182 mptval(icomp,2) = mptval(icomp,2) * debye 183 mptval(icomp,3) = mptval(icomp,1) + mptval(icomp,2) 184 enddo 185 dipol = dipol*debye 186 dipefc = dipefc*debye 187 diptot = diptot*debye 188 write(luout,2996) dipol,mptval(1,1),mptval(1,2),mptval(2,1), 189 & mptval(2,2),mptval(3,1),mptval(3,2),dipefc, 190 & diptot 191 write(luout,2994) 192 if (.not.rtdb_put(rtdb,'prop:dipval',mt_dbl,3,mptval(1,3))) 193 & call errquit('prop: rtdb_put of dipval failed',555, 194 & RTDB_ERR) 195c 196c Done, goto cleanup and sync with other nodes 197c 198 goto 999 199c 200c Quadrupole output 201c 202c Ordering in mptval (from defNxyz): 203c xx, xy, xz, yy, yz, zz 204c 205c First < r**2 > = diamagnetic susceptibility 206c 207 102 write(luout,3992) rsquar 208c 209c Second moments 210c 211 write(luout,3999) 212 write(luout,3996) mptval(1,1),mptval(1,2),mptval(1,1)+mptval(1,2), 213 1 mptval(4,1),mptval(4,2),mptval(4,1)+mptval(4,2), 214 2 mptval(6,1),mptval(6,2),mptval(6,1)+mptval(6,2), 215 3 mptval(2,1),mptval(2,2),mptval(2,1)+mptval(2,2), 216 4 mptval(3,1),mptval(3,2),mptval(3,1)+mptval(3,2), 217 5 mptval(5,1),mptval(5,2),mptval(5,1)+mptval(5,2) 218c 219 do j = 1, 2 220 do i = 1, 6 221 mptval2(i,j) = mptval(i,j)*buck 222 enddo 223 enddo 224c 225 write(luout,3995) 226 write(luout,3996) 227 1 mptval2(1,1),mptval2(1,2),mptval2(1,1)+mptval2(1,2), 228 2 mptval2(4,1),mptval2(4,2),mptval2(4,1)+mptval2(4,2), 229 3 mptval2(6,1),mptval2(6,2),mptval2(6,1)+mptval2(6,2), 230 4 mptval2(2,1),mptval2(2,2),mptval2(2,1)+mptval2(2,2), 231 5 mptval2(3,1),mptval2(3,2),mptval2(3,1)+mptval2(3,2), 232 6 mptval2(5,1),mptval2(5,2),mptval2(5,1)+mptval2(5,2) 233c 234c Quadrupole moments 235c 236 rsquar = mptval(1,1)+mptval(4,1)+mptval(6,1) 237 mptval(1,1)=(3.0d0*mptval(1,1)-rsquar)/2.0d0 238 mptval(4,1)=(3.0d0*mptval(4,1)-rsquar)/2.0d0 239 mptval(6,1)=(3.0d0*mptval(6,1)-rsquar)/2.0d0 240 mptval(2,1)= 3.0d0*mptval(2,1)/2.0d0 241 mptval(3,1)= 3.0d0*mptval(3,1)/2.0d0 242 mptval(5,1)= 3.0d0*mptval(5,1)/2.0d0 243 rsquar = mptval(1,2)+mptval(4,2)+mptval(6,2) 244 mptval(1,2)=(3.0d0*mptval(1,2)-rsquar)/2.0d0 245 mptval(4,2)=(3.0d0*mptval(4,2)-rsquar)/2.0d0 246 mptval(6,2)=(3.0d0*mptval(6,2)-rsquar)/2.0d0 247 mptval(2,2)= 3.0d0*mptval(2,2)/2.0d0 248 mptval(3,2)= 3.0d0*mptval(3,2)/2.0d0 249 mptval(5,2)= 3.0d0*mptval(5,2)/2.0d0 250 write(luout,3997) 251 write(luout,3996) mptval(1,1),mptval(1,2),mptval(1,1)+mptval(1,2), 252 1 mptval(4,1),mptval(4,2),mptval(4,1)+mptval(4,2), 253 2 mptval(6,1),mptval(6,2),mptval(6,1)+mptval(6,2), 254 3 mptval(2,1),mptval(2,2),mptval(2,1)+mptval(2,2), 255 4 mptval(3,1),mptval(3,2),mptval(3,1)+mptval(3,2), 256 5 mptval(5,1),mptval(5,2),mptval(5,1)+mptval(5,2) 257c 258 do j = 1, 2 259 do i = 1, 6 260 mptval2(i,j) = mptval(i,j)*buck 261 enddo 262 enddo 263c 264 write(luout,3994) 265 write(luout,3996) 266 1 mptval2(1,1),mptval2(1,2),mptval2(1,1)+mptval2(1,2), 267 2 mptval2(4,1),mptval2(4,2),mptval2(4,1)+mptval2(4,2), 268 3 mptval2(6,1),mptval2(6,2),mptval2(6,1)+mptval2(6,2), 269 4 mptval2(2,1),mptval2(2,2),mptval2(2,1)+mptval2(2,2), 270 5 mptval2(3,1),mptval2(3,2),mptval2(3,1)+mptval2(3,2), 271 6 mptval2(5,1),mptval2(5,2),mptval2(5,1)+mptval2(5,2) 272 write(luout,3993) 273c 274c Done, goto cleanup and sync with other nodes 275c 276 goto 999 277c 278c Octupole output 279c 280c Ordering in mptval (from defNxyz): 281c xxx, xxy, xxz, yyx, xyz, zzx, yyy, zzy, zzz 282c 283c Third moments 284c 285 103 write(luout,9999) 286 write(luout,9996) 287 1 mptval(1,1),mptval(1,2),mptval(1,1)+mptval(1,2), 288 2 mptval(7,1),mptval(7,2),mptval(7,1)+mptval(7,2), 289 3 mptval(10,1),mptval(10,2),mptval(10,1)+mptval(10,2), 290 4 mptval(2,1),mptval(2,2),mptval(2,1)+mptval(2,2), 291 5 mptval(3,1),mptval(3,2),mptval(3,1)+mptval(3,2), 292 6 mptval(4,1),mptval(4,2),mptval(4,1)+mptval(4,2), 293 7 mptval(8,1),mptval(8,2),mptval(8,1)+mptval(8,2), 294 8 mptval(6,1),mptval(6,2),mptval(6,1)+mptval(6,2), 295 9 mptval(9,1),mptval(9,2),mptval(9,1)+mptval(9,2), 296 1 mptval(5,1),mptval(5,2),mptval(5,1)+mptval(5,2) 297c 298 do j = 1, 2 299 do i = 1, 10 300 mptval2(i,j) = mptval(i,j)*octo 301 enddo 302 enddo 303c 304 write(luout,9995) 305 write(luout,9996) 306 1 mptval2(1,1),mptval2(1,2),mptval2(1,1)+mptval2(1,2), 307 2 mptval2(7,1),mptval2(7,2),mptval2(7,1)+mptval2(7,2), 308 3 mptval2(10,1),mptval2(10,2),mptval2(10,1)+mptval2(10,2), 309 4 mptval2(2,1),mptval2(2,2),mptval2(2,1)+mptval2(2,2), 310 5 mptval2(3,1),mptval2(3,2),mptval2(3,1)+mptval2(3,2), 311 6 mptval2(4,1),mptval2(4,2),mptval2(4,1)+mptval2(4,2), 312 7 mptval2(8,1),mptval2(8,2),mptval2(8,1)+mptval2(8,2), 313 8 mptval2(6,1),mptval2(6,2),mptval2(6,1)+mptval2(6,2), 314 9 mptval2(9,1),mptval2(9,2),mptval2(9,1)+mptval2(9,2), 315 1 mptval2(5,1),mptval2(5,2),mptval2(5,1)+mptval2(5,2) 316c 317c Octupole moments 318c 319 do i=1,2 320 dum1=mptval(1,i) 321 dum2=mptval(4,i) 322 mptval(1,i)=mptval(1,i)-3.0d0*(mptval(4,i)+mptval(6,i))/2.0d0 323 mptval(4,i)=(4.0d0*mptval(4,i)-dum1-mptval(6,i))/2.0d0 324 mptval(6,i)=(4.0d0*mptval(6,i)-dum1-dum2)/2.0d0 325 dum1=mptval(7,i) 326 dum2=mptval(2,i) 327 mptval(7,i)=mptval(7,i)-3.0d0*(mptval(2,i)+mptval(9,i))/2.0d0 328 mptval(2,i)=(4.0d0*mptval(2,i)-dum1-mptval(9,i))/2.0d0 329 mptval(9,i)=(4.0d0*mptval(9,i)-dum1-dum2)/2.0d0 330 dum1=mptval(10,i) 331 dum2=mptval(3,i) 332 mptval(10,i)=mptval(10,i)-3.0d0*(mptval(3,i)+mptval(8,i))/2.0d0 333 mptval(3,i)=(4.0d0*mptval(3,i)-dum1-mptval(8,i))/2.0d0 334 mptval(8,i)=(4.0d0*mptval(8,i)-dum1-dum2)/2.0d0 335 mptval(5,i)=5.0d0*mptval(5,i)/2.0d0 336 enddo 337 write(luout,9997) 338 write(luout,9996) 339 1 mptval(1,1),mptval(1,2),mptval(1,1)+mptval(1,2), 340 2 mptval(7,1),mptval(7,2),mptval(7,1)+mptval(7,2), 341 3 mptval(10,1),mptval(10,2),mptval(10,1)+mptval(10,2), 342 4 mptval(2,1),mptval(2,2),mptval(2,1)+mptval(2,2), 343 5 mptval(3,1),mptval(3,2),mptval(3,1)+mptval(3,2), 344 6 mptval(4,1),mptval(4,2),mptval(4,1)+mptval(4,2), 345 7 mptval(8,1),mptval(8,2),mptval(8,1)+mptval(8,2), 346 8 mptval(6,1),mptval(6,2),mptval(6,1)+mptval(6,2), 347 9 mptval(9,1),mptval(9,2),mptval(9,1)+mptval(9,2), 348 1 mptval(5,1),mptval(5,2),mptval(5,1)+mptval(5,2) 349c 350 do j = 1, 2 351 do i = 1, 10 352 mptval2(i,j) = mptval(i,j)*octo 353 enddo 354 enddo 355c 356 write(luout,9994) 357 write(luout,9996) 358 1 mptval2(1,1),mptval2(1,2),mptval2(1,1)+mptval2(1,2), 359 2 mptval2(7,1),mptval2(7,2),mptval2(7,1)+mptval2(7,2), 360 3 mptval2(10,1),mptval2(10,2),mptval2(10,1)+mptval2(10,2), 361 4 mptval2(2,1),mptval2(2,2),mptval2(2,1)+mptval2(2,2), 362 5 mptval2(3,1),mptval2(3,2),mptval2(3,1)+mptval2(3,2), 363 6 mptval2(4,1),mptval2(4,2),mptval2(4,1)+mptval2(4,2), 364 7 mptval2(8,1),mptval2(8,2),mptval2(8,1)+mptval2(8,2), 365 8 mptval2(6,1),mptval2(6,2),mptval2(6,1)+mptval2(6,2), 366 9 mptval2(9,1),mptval2(9,2),mptval2(9,1)+mptval2(9,2), 367 1 mptval2(5,1),mptval2(5,2),mptval2(5,1)+mptval2(5,2) 368 write(luout,9993) 369c 370c ----- release MA memory blocks ----- 371c 372 999 if (.not.ma_pop_stack(l_zanpt)) call errquit 373 & ('hnd_mtpole, ma_pop_stack of l_zanpt failed',911,MA_ERR) 374 if (.not.ma_pop_stack(l_xyzpt)) call errquit 375 & ('hnd_mtpole, ma_pop_stack of l_xyzpt failed',911,MA_ERR) 376c 377c Synchronize all nodes (as only node 0 was writing) 378c Reset rtdb access to parallel 379c 380 1000 call ga_sync() 381 status = rtdb_parallel(.true.) 382c 383 do i = 1, ndens 384 if (.not.ga_destroy(g_dens(i))) call 385 & errquit('mtpole: ga_destroy failed g_dens',0,GA_ERR) 386 enddo 387c 388c Terminate integrals 389c 390 call schwarz_tidy() 391 call int_terminate() 392c 393 return 394c 395c Formatting dipole 396c 397 2001 format(/,10x,13(1h-),/,10x,'Dipole Moment',/,10x,13(1h-)) 398 2996 format(/,2x,' Dipole moment',f20.10,' Debye(s)', 399 2 /,12x,' DMX',f20.10,' DMXEFC',f20.10, 400 3 /,12x,' DMY',f20.10,' DMYEFC',f20.10, 401 4 /,12x,' DMZ',f20.10,' DMZEFC',f20.10, 402 5 /,2x,' -EFC- dipole ',f20.10,' DEBYE(S)', 403 6 /,2x,' Total dipole ',f20.10,' DEBYE(S)') 404 2995 format(/,2x,' Dipole moment',f20.10,' A.U.', 405 2 /,12x,' DMX',f20.10,' DMXEFC',f20.10, 406 3 /,12x,' DMY',f20.10,' DMYEFC',f20.10, 407 4 /,12x,' DMZ',f20.10,' DMZEFC',f20.10, 408 5 /, 2x,' -EFC- dipole ',f20.10,' A.U.', 409 6 /, 2x,' Total dipole ',f20.10,' A.U.') 410 2994 format(/,' 1 a.u. = 2.541766 Debyes ') 411c 412c Formatting quadrupole 413c 414 3001 format(/,10x,17(1h-),/,10x,'Quadrupole Moment',/,10x,17(1h-)) 415 3999 format(/,2x,' Second moments in atomic units') 416 3997 format(/,2x,' Quadrupole moments in atomic units') 417 3995 format(/,2x,' Second moments in buckingham(s)') 418 3994 format(/,2x,' Quadrupole moments in buckingham(s)') 419 3996 format(/,2x,' Component',' Electronic+nuclear',4x, 420 2 ' Point charges',12x,' Total',/,2x,74(1h-), 421 2 /,2x,' XX ',f20.10,2x,f20.10,2x,f20.10, 422 2 /,2x,' YY ',f20.10,2x,f20.10,2x,f20.10, 423 2 /,2x,' ZZ ',f20.10,2x,f20.10,2x,f20.10, 424 2 /,2x,' XY ',f20.10,2x,f20.10,2x,f20.10, 425 2 /,2x,' XZ ',f20.10,2x,f20.10,2x,f20.10, 426 2 /,2x,' YZ ',f20.10,2x,f20.10,2x,f20.10) 427 3993 format(/,' 1 a.u. = 1.344911 Buckinghams ', 428 1 '= 1.344911 10**(-26) esu*cm**2 ') 429 3992 format(/,' < R**2 > = ',f10.6,' a.u. ', 430 1 ' ( 1 a.u. = 0.280023 10**(-16) cm**2 ) ',/, 431 2 ' ( also called diamagnetic susceptibility ) ') 432c 433c Formatted printing octupole moments 434c 435 9001 format(/,10x,15(1h-),/,10x,'Octupole Moment',/,10x,15(1h-)) 436 9999 format(/,2x,' Third moments in atomic units') 437 9997 format(/,2x,' Octupole moments in atomic units') 438 9995 format(/,2x,' Third moments in 10**(-34) esu*cm**3') 439 9994 format(/,2x,' Octupole moments in 10**(-34) esu*cm**3') 440 9996 format(/,2x,' Component',' Electronic+nuclear',4x, 441 2 ' Point charges',12x,' Total',/,2x,74(1h-), 442 2 /,2x,' XXX ',f20.10,2x,f20.10,2x,f20.10, 443 2 /,2x,' YYY ',f20.10,2x,f20.10,2x,f20.10, 444 2 /,2x,' ZZZ ',f20.10,2x,f20.10,2x,f20.10, 445 2 /,2x,' XXY ',f20.10,2x,f20.10,2x,f20.10, 446 2 /,2x,' XXZ ',f20.10,2x,f20.10,2x,f20.10, 447 2 /,2x,' YYX ',f20.10,2x,f20.10,2x,f20.10, 448 2 /,2x,' YYZ ',f20.10,2x,f20.10,2x,f20.10, 449 2 /,2x,' ZZX ',f20.10,2x,f20.10,2x,f20.10, 450 2 /,2x,' ZZY ',f20.10,2x,f20.10,2x,f20.10, 451 2 /,2x,' XYZ ',f20.10,2x,f20.10,2x,f20.10) 452 9993 format(/,' 1 a.u. = 0.711688 10**(-34) esu*cm**3 ') 453 end 454c 455 subroutine prp_cent(rtdb,geom,cm) 456c 457c Routine calculates coordinates of the expansion center 458c depending on what kind of expansion the user chooses 459c 1: center of charge (default) 460c 2: center of mass 461c 3: origin 462c 4: arbitrary point 463c 464 implicit none 465#include "errquit.fh" 466#include "geom.fh" 467#include "rtdb.fh" 468#include "global.fh" 469#include "mafdecls.fh" 470#include "stdio.fh" 471c 472 integer rtdb ! [input] rtdb handle 473 integer geom ! [input] geometry handle 474 double precision cm(3) ! [output] expansion center 475c 476 integer i, i_expan 477 double precision scale 478c 479c Read in point of expansion 480c 481 if (.not. rtdb_get(rtdb, 'prop:center', mt_int, 1, i_expan)) 482 $ i_expan = 1 ! default of center of charge 483c 484 if (i_expan.eq.4) then ! arbitrary point of expansion 485 if (.not. rtdb_get(rtdb, 'prop:center_val',mt_dbl,3,cm)) 486 $ call errquit('prp_cent: rtdb_get center_val failed',555, 487 & RTDB_ERR) 488c 489c get scale factor based on geometry UNITS keyword 490c 491 if (.not. geom_get_user_scale(geom,scale)) 492 $ call errquit('prp_cent: trouble getting units scale',555, 493 & GEOM_ERR) 494 do i = 1, 3 ! scale according to geometry units 495 cm(i) = cm(i)*scale 496 end do 497 if(ga_nodeid().eq.0) write(luout,5994) cm(1),cm(2),cm(3) 498c 499 elseif (i_expan.eq.3) then ! origin 500 do i = 1, 3 501 cm(i) = 0.d0 502 end do 503 if (ga_nodeid().eq.0) write(luout,5995) cm(1),cm(2),cm(3) 504c 505 elseif (i_expan.eq.2) then ! center of mass 506 if (.not.geom_center_of_mass(geom,cm)) call 507 & errquit('prp_cent: geom_center_of_mass failed',0,GEOM_ERR) 508 if (ga_nodeid().eq.0) write(luout,5996) cm(1),cm(2),cm(3) 509c 510 elseif (i_expan.eq.1) then ! center of charge 511 if (.not.geom_center_of_charge(geom,cm)) call 512 & errquit('prp_cent: geom_center_of_charge failed',0,GEOM_ERR) 513 if (ga_nodeid().eq.0) write(luout,5997) cm(1),cm(2),cm(3) 514 else 515 call errquit('prp_cent: improper center value',555, INPUT_ERR) 516 endif 517c 518 return 519 5994 format(/,' Center of expansion (in au) is the arbitrary point', 520 1 /,8x,' X = ',f15.7,' Y = ',f15.7,' Z = ',f15.7) 521 5995 format(/,' Center of expansion (in au) is the origin', 522 1 /,8x,' X = ',f15.7,' Y = ',f15.7,' Z = ',f15.7) 523 5996 format(/,' Center of mass (in au) is the expansion point', 524 1 /,8x,' X = ',f15.7,' Y = ',f15.7,' Z = ',f15.7) 525 5997 format(/,' Center of charge (in au) is the expansion point', 526 1 /,8x,' X = ',f15.7,' Y = ',f15.7,' Z = ',f15.7) 527 end 528c 529c 530c 531c 532c 533c Jeff wrote a new version which returns mtpval - in atomic units - and 534c prints less stuff, for use in the TCE. 535c 536 subroutine hnd_mtpole2(rtdb,basis,geom,mpole,mtpval) 537c 538c $Id$ 539c 540c This routine calculates the multipole moments around 541c choosen center of expansion 542c 543c mpole = 1 : Dipole 544c mpole = 2 : Quadrupole 545c mpole = 3 : Octupole 546c 547 implicit none 548#include "errquit.fh" 549#include "global.fh" 550#include "mafdecls.fh" 551#include "nwc_const.fh" 552#include "stdio.fh" 553#include "geom.fh" 554#include "rtdb.fh" 555#include "cosmo.fh" 556c 557 integer rtdb ! [input] rtdb handle 558 integer basis ! [input] basis handle 559 integer geom ! [input] geometry handle 560 integer mpole ! [input] order of multipole 561 integer max_pole, maxcomp 562 parameter (max_pole = 3, maxcomp = (max_pole+1)*(max_pole+2)/2) 563c 564 double precision cm(3), mtpval(maxcomp,3) 565 double precision cx, cy, cz 566 double precision dipol, dipefc, diptot, rsquar, dum1, dum2 567 logical status 568 integer Nxyz(3) 569 integer ncomp, nat, iat, icomp, i 570 integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt 571 integer nefc, l_efcc, k_efcc, l_efcz, k_efcz 572 integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2) 573 character*3 scftyp 574 character*16 at_tag 575 if (ga_nodeid().eq.0.and.mpole.eq.1) write(luout,2001) 576 if (ga_nodeid().eq.0.and.mpole.eq.2) write(luout,3001) 577 if (ga_nodeid().eq.0.and.mpole.eq.3) write(luout,9001) 578c 579c Get density matrix 580c 581 call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp, 582 & nclosed,nopen,nvirt) 583c 584c Determine number of components and zero mtpval 585c 586 ncomp = (mpole+1)*(mpole+2)/2 587c 588c It would appear dcopy screws things up be zeroing the wrong area. 589c This bothers me but I do not care now that it is fixed. 590c call dcopy(2*maxcomp,0.0d0,0,mtpval,1) 591c 592 do icomp = 1, ncomp 593 mtpval(icomp,1) = 0.0d0 594 mtpval(icomp,2) = 0.0d0 595 mtpval(icomp,3) = 0.0d0 596 enddo 597c 598c Get center of expansion 599c 600 call prp_cent(rtdb,geom,cm) 601c 602c Electronic contribution 603c 604 call hnd_mtpcon(basis,geom,g_dens(ndens),mtpval,mpole,cm) 605c 606c Only node 0 writing to rtdb 607c 608 status = rtdb_parallel(.false.) 609c 610 if (ga_nodeid().gt.0) goto 1000 611c 612c Allocate some local memory 613c 614 if (.not.geom_ncent(geom,nat)) call 615 & errquit('hnd_efgmap: geom_ncent failed',911,GEOM_ERR) 616 if (.not. ma_push_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 617 & call errquit('hnd_mtpole2: ma failed',911,MA_ERR) 618 if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 619 & call errquit('hnd_mtpole2: ma failed',911,MA_ERR) 620c 621 do iat=1,nat 622 if (.not.geom_cent_get(geom,iat,at_tag, 623 & dbl_mb(k_xyzpt+3*(iat-1)),dbl_mb(k_zanpt+iat-1))) call 624 & errquit('hnd_mtpole2: geom_cent_get failed',911,GEOM_ERR) 625 enddo ! iat 626c 627c For quadrupole case, need electronic contribution for 628c diamagnetic susceptibility 629c 630 if (mpole.eq.2) rsquar=mtpval(1,1)+mtpval(4,1)+mtpval(6,1) 631c 632c Nuclear contribution 633c 634 do icomp = 1, ncomp 635 mtpval(icomp,1) = -mtpval(icomp,1) 636 call getNxyz(mpole,icomp,Nxyz) 637 do iat = 1, nat 638 cx = 1d0 639 if(Nxyz(1).ne.0) 640 & cx = (dbl_mb(k_xyzpt +3*(iat-1))-cm(1))**Nxyz(1) 641 cy = 1d0 642 if(Nxyz(2).ne.0) 643 & cy = (dbl_mb(k_xyzpt+1+3*(iat-1))-cm(2))**Nxyz(2) 644 cz = 1d0 645 if(Nxyz(3).ne.0) 646 & cz = (dbl_mb(k_xyzpt+2+3*(iat-1))-cm(3))**Nxyz(3) 647 mtpval(icomp,1) = mtpval(icomp,1) + 648 & dbl_mb(k_zanpt+iat-1) * cx * cy * cz 649 enddo ! iat 650 enddo ! icomp 651c 652c Print output for 653c mpole = 1 : Dipole 654c mpole = 2 : Quadrupole 655c mpole = 3 : Octupole 656c 657 goto (101,102,103) mpole 658c 659c Dipole output 660c 661 101 continue 662c 663 do icomp = 1, ncomp 664 mtpval(icomp,1) = -1.0*mtpval(icomp,1) 665 enddo 666c 667 goto 999 668c 669c Quadrupole output 670c 671c Ordering in mtpval (from defNxyz): 672c xx, xy, xz, yy, yz, zz 673c 674 102 continue 675c 676c Quadrupole moments 677c 678 rsquar = mtpval(1,1)+mtpval(4,1)+mtpval(6,1) 679 mtpval(1,1)=(3.0d0*mtpval(1,1)-rsquar)/2.0d0 680 mtpval(4,1)=(3.0d0*mtpval(4,1)-rsquar)/2.0d0 681 mtpval(6,1)=(3.0d0*mtpval(6,1)-rsquar)/2.0d0 682 mtpval(2,1)= 3.0d0*mtpval(2,1)/2.0d0 683 mtpval(3,1)= 3.0d0*mtpval(3,1)/2.0d0 684 mtpval(5,1)= 3.0d0*mtpval(5,1)/2.0d0 685 rsquar = mtpval(1,2)+mtpval(4,2)+mtpval(6,2) 686 mtpval(1,2)=(3.0d0*mtpval(1,2)-rsquar)/2.0d0 687 mtpval(4,2)=(3.0d0*mtpval(4,2)-rsquar)/2.0d0 688 mtpval(6,2)=(3.0d0*mtpval(6,2)-rsquar)/2.0d0 689 mtpval(2,2)= 3.0d0*mtpval(2,2)/2.0d0 690 mtpval(3,2)= 3.0d0*mtpval(3,2)/2.0d0 691 mtpval(5,2)= 3.0d0*mtpval(5,2)/2.0d0 692c 693 goto 999 694c 695c Octupole output 696c 697c Ordering in mtpval (from defNxyz): 698c xxx, xxy, xxz, yyx, xyz, zzx, yyy, zzy, zzz 699c 700 103 continue 701c 702c Octupole moments 703c 704 do i=1,2 705 dum1=mtpval(1,i) 706 dum2=mtpval(4,i) 707 mtpval(1,i)=mtpval(1,i)-3.0d0*(mtpval(4,i)+mtpval(6,i))/2.0d0 708 mtpval(4,i)=(4.0d0*mtpval(4,i)-dum1-mtpval(6,i))/2.0d0 709 mtpval(6,i)=(4.0d0*mtpval(6,1)-dum1-dum2)/2.0d0 710 dum1=mtpval(7,i) 711 dum2=mtpval(2,i) 712 mtpval(7,i)=mtpval(7,i)-3.0d0*(mtpval(2,i)+mtpval(9,i))/2.0d0 713 mtpval(2,i)=(4.0d0*mtpval(2,i)-dum1-mtpval(9,i))/2.0d0 714 mtpval(9,i)=(4.0d0*mtpval(9,i)-dum1-dum2)/2.0d0 715 dum1=mtpval(10,i) 716 dum2=mtpval(3,i) 717 mtpval(10,i)=mtpval(10,i)-3.0d0*(mtpval(3,i)+mtpval(8,i))/2.0d0 718 mtpval(3,i)=(4.0d0*mtpval(3,i)-dum1-mtpval(8,i))/2.0d0 719 mtpval(8,i)=(4.0d0*mtpval(8,i)-dum1-dum2)/2.0d0 720 mtpval(5,i)=5.0d0*mtpval(5,i)/2.0d0 721 enddo 722c 723c ----- release MA memory blocks ----- 724c 725 999 if (.not.ma_pop_stack(l_zanpt)) call errquit 726 & ('hnd_mtpole2, ma_pop_stack of l_zanpt failed',911,MA_ERR) 727 if (.not.ma_pop_stack(l_xyzpt)) call errquit 728 & ('hnd_mtpole2, ma_pop_stack of l_xyzpt failed',911,MA_ERR) 729c 730c Synchronize all nodes (as only node 0 was writing) 731c Reset rtdb access to parallel 732c 733 1000 call ga_sync() 734 status = rtdb_parallel(.true.) 735c 736 do i = 1, ndens 737 if (.not.ga_destroy(g_dens(i))) call 738 & errquit('mtpole: ga_destroy failed g_dens',0,GA_ERR) 739 enddo 740c 741 return 742c 743 2001 format(/,10x,13(1h-),/,10x,'Dipole Moment',/,10x,13(1h-)) 744 3001 format(/,10x,17(1h-),/,10x,'Quadrupole Moment',/,10x,17(1h-)) 745 9001 format(/,10x,15(1h-),/,10x,'Octupole Moment',/,10x,15(1h-)) 746 end 747 748 subroutine hnd_mtpole_nuclear(rtdb,geom,mpole,mtpval) 749c 750c $Id$ 751c 752c This routine calculates the multipole moments around 753c choosen center of expansion 754c 755c mpole = 1 : Dipole 756c mpole = 2 : Quadrupole 757c mpole = 3 : Octupole 758c 759 implicit none 760#include "errquit.fh" 761#include "global.fh" 762#include "mafdecls.fh" 763#include "nwc_const.fh" 764#include "stdio.fh" 765#include "geom.fh" 766#include "rtdb.fh" 767#include "cosmo.fh" 768c 769 integer rtdb ! [input] rtdb handle 770 integer geom ! [input] geometry handle 771 integer mpole ! [input] order of multipole 772 integer max_pole, maxcomp 773 parameter (max_pole = 3, maxcomp = (max_pole+1)*(max_pole+2)/2) 774c 775 double precision cm(3), mtpval(maxcomp,3) 776 double precision cx, cy, cz 777 double precision dipol, dipefc, diptot, rsquar, dum1, dum2 778 logical status 779 integer Nxyz(3) 780 integer ncomp, nat, iat, icomp, i 781 integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt 782 integer nefc, l_efcc, k_efcc, l_efcz, k_efcz 783 integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2) 784 character*3 scftyp 785 character*16 at_tag 786 if (ga_nodeid().eq.0.and.mpole.eq.1) write(luout,2001) 787 if (ga_nodeid().eq.0.and.mpole.eq.2) write(luout,3001) 788 if (ga_nodeid().eq.0.and.mpole.eq.3) write(luout,9001) 789c 790c Determine number of components and zero mtpval 791c 792 ncomp = (mpole+1)*(mpole+2)/2 793c 794c It would appear dcopy screws things up be zeroing the wrong area. 795c This bothers me but I do not care now that it is fixed. 796c call dcopy(2*maxcomp,0.0d0,0,mtpval,1) 797c 798 do icomp = 1, ncomp 799 mtpval(icomp,1) = 0.0d0 800 mtpval(icomp,2) = 0.0d0 801 mtpval(icomp,3) = 0.0d0 802 enddo 803c 804c Get center of expansion 805c 806 call prp_cent(rtdb,geom,cm) 807c 808c Only node 0 writing to rtdb 809c 810 status = rtdb_parallel(.false.) 811c 812 if (ga_nodeid().gt.0) goto 1000 813c 814c Allocate some local memory 815c 816 if (.not.geom_ncent(geom,nat)) call 817 & errquit('hnd_efgmap: geom_ncent failed',911,GEOM_ERR) 818 if (.not. ma_push_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 819 & call errquit('hnd_mtpole2: ma failed',911,MA_ERR) 820 if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 821 & call errquit('hnd_mtpole2: ma failed',911,MA_ERR) 822c 823 do iat=1,nat 824 if (.not.geom_cent_get(geom,iat,at_tag, 825 & dbl_mb(k_xyzpt+3*(iat-1)),dbl_mb(k_zanpt+iat-1))) call 826 & errquit('hnd_mtpole2: geom_cent_get failed',911,GEOM_ERR) 827 enddo ! iat 828c 829c For quadrupole case, need electronic contribution for 830c diamagnetic susceptibility 831c 832 if (mpole.eq.2) rsquar=mtpval(1,1)+mtpval(4,1)+mtpval(6,1) 833c 834c Nuclear contribution 835c 836 do icomp = 1, ncomp 837 mtpval(icomp,1) = -mtpval(icomp,1) 838 call getNxyz(mpole,icomp,Nxyz) 839 do iat = 1, nat 840 cx = 1d0 841 if(Nxyz(1).ne.0) 842 & cx = (dbl_mb(k_xyzpt +3*(iat-1))-cm(1))**Nxyz(1) 843 cy = 1d0 844 if(Nxyz(2).ne.0) 845 & cy = (dbl_mb(k_xyzpt+1+3*(iat-1))-cm(2))**Nxyz(2) 846 cz = 1d0 847 if(Nxyz(3).ne.0) 848 & cz = (dbl_mb(k_xyzpt+2+3*(iat-1))-cm(3))**Nxyz(3) 849 mtpval(icomp,1) = mtpval(icomp,1) + 850 & dbl_mb(k_zanpt+iat-1) * cx * cy * cz 851 enddo ! iat 852 enddo ! icomp 853c 854c Print output for 855c mpole = 1 : Dipole 856c mpole = 2 : Quadrupole 857c mpole = 3 : Octupole 858c 859 goto (101,102,103) mpole 860c 861c Dipole output 862c 863 101 continue 864c 865 do icomp = 1, ncomp 866 mtpval(icomp,1) = -1.0*mtpval(icomp,1) 867 enddo 868c 869 goto 999 870c 871c Quadrupole output 872c 873c Ordering in mtpval (from defNxyz): 874c xx, xy, xz, yy, yz, zz 875c 876 102 continue 877c 878c Quadrupole moments 879c 880 rsquar = mtpval(1,1)+mtpval(4,1)+mtpval(6,1) 881 mtpval(1,1)=(3.0d0*mtpval(1,1)-rsquar)/2.0d0 882 mtpval(4,1)=(3.0d0*mtpval(4,1)-rsquar)/2.0d0 883 mtpval(6,1)=(3.0d0*mtpval(6,1)-rsquar)/2.0d0 884 mtpval(2,1)= 3.0d0*mtpval(2,1)/2.0d0 885 mtpval(3,1)= 3.0d0*mtpval(3,1)/2.0d0 886 mtpval(5,1)= 3.0d0*mtpval(5,1)/2.0d0 887 rsquar = mtpval(1,2)+mtpval(4,2)+mtpval(6,2) 888 mtpval(1,2)=(3.0d0*mtpval(1,2)-rsquar)/2.0d0 889 mtpval(4,2)=(3.0d0*mtpval(4,2)-rsquar)/2.0d0 890 mtpval(6,2)=(3.0d0*mtpval(6,2)-rsquar)/2.0d0 891 mtpval(2,2)= 3.0d0*mtpval(2,2)/2.0d0 892 mtpval(3,2)= 3.0d0*mtpval(3,2)/2.0d0 893 mtpval(5,2)= 3.0d0*mtpval(5,2)/2.0d0 894c 895 goto 999 896c 897c Octupole output 898c 899c Ordering in mtpval (from defNxyz): 900c xxx, xxy, xxz, yyx, xyz, zzx, yyy, zzy, zzz 901c 902 103 continue 903c 904c Octupole moments 905c 906 do i=1,2 907 dum1=mtpval(1,i) 908 dum2=mtpval(4,i) 909 mtpval(1,i)=mtpval(1,i)-3.0d0*(mtpval(4,i)+mtpval(6,i))/2.0d0 910 mtpval(4,i)=(4.0d0*mtpval(4,i)-dum1-mtpval(6,i))/2.0d0 911 mtpval(6,i)=(4.0d0*mtpval(6,1)-dum1-dum2)/2.0d0 912 dum1=mtpval(7,i) 913 dum2=mtpval(2,i) 914 mtpval(7,i)=mtpval(7,i)-3.0d0*(mtpval(2,i)+mtpval(9,i))/2.0d0 915 mtpval(2,i)=(4.0d0*mtpval(2,i)-dum1-mtpval(9,i))/2.0d0 916 mtpval(9,i)=(4.0d0*mtpval(9,i)-dum1-dum2)/2.0d0 917 dum1=mtpval(10,i) 918 dum2=mtpval(3,i) 919 mtpval(10,i)=mtpval(10,i)-3.0d0*(mtpval(3,i)+mtpval(8,i))/2.0d0 920 mtpval(3,i)=(4.0d0*mtpval(3,i)-dum1-mtpval(8,i))/2.0d0 921 mtpval(8,i)=(4.0d0*mtpval(8,i)-dum1-dum2)/2.0d0 922 mtpval(5,i)=5.0d0*mtpval(5,i)/2.0d0 923 enddo 924c 925c ----- release MA memory blocks ----- 926c 927 999 if (.not.ma_pop_stack(l_zanpt)) call errquit 928 & ('hnd_mtpole2, ma_pop_stack of l_zanpt failed',911,MA_ERR) 929 if (.not.ma_pop_stack(l_xyzpt)) call errquit 930 & ('hnd_mtpole2, ma_pop_stack of l_xyzpt failed',911,MA_ERR) 931c 932c Synchronize all nodes (as only node 0 was writing) 933c Reset rtdb access to parallel 934c 935 1000 call ga_sync() 936 status = rtdb_parallel(.true.) 937c 938 do i = 1, ndens 939 if (.not.ga_destroy(g_dens(i))) call 940 & errquit('mtpole: ga_destroy failed g_dens',0,GA_ERR) 941 enddo 942c 943 return 944c 945 2001 format(/,10x,13(1h-),/,10x,'Dipole Moment',/,10x,13(1h-)) 946 3001 format(/,10x,17(1h-),/,10x,'Quadrupole Moment',/,10x,17(1h-)) 947 9001 format(/,10x,15(1h-),/,10x,'Octupole Moment',/,10x,15(1h-)) 948 end 949