1* 2* $Id$ 3* 4 5* *********************************************** 6* * * 7* * pspw_dplot * 8* * * 9* *********************************************** 10 11 logical function pspw_dplot(rtdb) 12 implicit none 13 integer rtdb 14 15#include "global.fh" 16#include "bafdecls.fh" 17#include "btdb.fh" 18#include "errquit.fh" 19 20 21 logical value,fractional 22 integer taskid,np 23 integer MASTER 24 parameter (MASTER=0) 25 26* **** local variables *** 27 28 integer count,smearoccupation 29 integer ngrid(3) 30 integer nfft3d,n2ft3d,npack1,mapping,mapping1d 31 integer ispin,ne(2),nemax 32 integer psi2(2),dn(2),psir(2),occ2(2),ncell(3) 33 real*8 cpu1,cpu2 34 real*8 position_tolerance,origin(3) 35 character*8 tag 36 integer psir2(2) 37 38* **** external functions **** 39 logical control_read,ion_init,control_balance 40 integer pack_nwave_all,control_mapping,control_np_orbital 41 integer control_ngrid,pack_nwave,control_version 42 integer control_mapping1d 43 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 44 real*8 lattice_unitg 45 character*12 control_boundry 46 47 external control_read,ion_init,control_balance 48 external pack_nwave_all,control_mapping,control_np_orbital 49 external control_ngrid,pack_nwave,control_version 50 external control_mapping1d 51 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 52 external lattice_unitg 53 external control_boundry 54 55* |************| 56*****************************| PROLOGUE |**************************** 57* |************| 58 59 value = .true. 60 61* **** get parallel variables **** 62 call Parallel_Init() 63 call Parallel_np(np) 64 call Parallel_taskid(taskid) 65 if (taskid.eq.MASTER) call current_second(cpu1) 66 67* ***** print out header **** 68 if (taskid.eq.MASTER) then 69 write(6,1000) 70 write(6,1010) 71 write(6,1020) 72 write(6,1010) 73 write(6,1030) 74 write(6,1010) 75 write(6,1035) 76 write(6,1010) 77 write(6,1040) 78 write(6,1010) 79 write(6,1000) 80 write(6,*) 81 call nwpw_message(1) 82 write(6,1110) 83 call util_flush(6) 84 end if 85 86* **** get position_tolerance **** 87 if (.not. 88 > btdb_get(rtdb,'pspw_dplot:position_tolerance', 89 > mt_dbl,1,position_tolerance)) 90 > position_tolerance = 0.0d0 91 92* **** get ncell **** 93 if (.not.btdb_get(rtdb,'pspw_dplot:ncell',mt_int,3,ncell)) then 94 ncell(1) = 0 95 ncell(2) = 0 96 ncell(3) = 0 97 end if 98 99* **** get origin **** 100 if (.not. 101 > btdb_get(rtdb,'pspw_dplot:origin', 102 > mt_dbl,3,origin)) then 103 origin(1) = 0.0d0 104 origin(2) = 0.0d0 105 origin(3) = 0.0d0 106 end if 107 108* **** read control file **** 109 value = control_read(4,rtdb) 110 ngrid(1) = control_ngrid(1) 111 ngrid(2) = control_ngrid(2) 112 ngrid(3) = control_ngrid(3) 113 mapping = control_mapping() 114 115 116* **** initialize D3dB data structure **** 117 call D3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 118 call D3dB_nfft3d(1,nfft3d) 119 n2ft3d = 2*nfft3d 120 121* ***** Initialize double D3dB data structure **** 122 if (control_version().eq.4) 123 > call D3dB_Init(2,2*ngrid(1),2*ngrid(2),2*ngrid(3),mapping) 124 125* **** initialize lattice and packing data structure **** 126 call lattice_init() 127 call G_init() 128 call mask_init() 129 call Pack_init() 130 call Pack_npack(1,npack1) 131 132* **** initialize ion data structure **** 133 value = ion_init(rtdb) 134 135* ***** allocate psi2 wavefunctions **** 136 call psi_get_ne_occupation(ispin,ne,smearoccupation) 137 if (smearoccupation.gt.0) then 138 fractional = .true. 139 else 140 fractional = .false. 141 end if 142 nemax = ne(1)+ne(2) 143 value = value.and. 144 > BA_alloc_get(mt_dcpl,npack1*(ne(1)+ne(2)), 145 > 'psi2',psi2(2),psi2(1)) 146 if (fractional) then 147 value = value.and. 148 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ2',occ2(2),occ2(1)) 149 end if 150 151 if (.not. value) call errquit('out of heap memory',0, MA_ERR) 152 153 mapping1d = control_mapping1d() 154 call Dne_init(ispin,ne,mapping1d) 155 156* ***** read psi2 wavefunctions **** 157 call psi_read(ispin,ne,dcpl_mb(psi2(1)), 158 > smearoccupation,dbl_mb(occ2(1))) 159 160* **** allocate other variables ***** 161 value = value.and. 162 > BA_alloc_get(mt_dbl,(4*nfft3d), 163 > 'dn',dn(2),dn(1)) 164 value = value.and. 165 > BA_alloc_get(mt_dcpl,nfft3d*(ne(1)+ne(2)), 166 > 'psir',psir(2),psir(1)) 167 value = value.and. 168 > BA_alloc_get(mt_dcpl,nfft3d*(ne(1)+ne(2)), 169 > 'psir2',psir2(2),psir2(1)) 170 171 if (.not. value) call errquit('pspw_dplot:out of heap memory',0, 172 & MA_ERR) 173 174 call ke_init() 175 if (control_version().eq.3) call coulomb_init() 176 if (control_version().eq.4) call coulomb2_init() 177 call strfac_init() 178 call phafac() 179 180* |**************************| 181****************** summary of input data ********************** 182* |**************************| 183 184 if (taskid.eq.MASTER) then 185 write(6,1111) np 186 if (mapping.eq.1) write(6,1112) 187 if (mapping.eq.2) write(6,1113) 188 if (mapping.eq.3) write(6,1118) 189 if (control_balance()) then 190 write(6,1114) 191 else 192 write(6,1116) 193 end if 194 write(6,1115) 195 write(6,1121) control_boundry(),control_version() 196 write(6,1220) ne(1),ne(ispin),' ( Fourier space)' 197 write(6,1224) ncell 198 write(6,1225) position_tolerance 199 write(6,1226) origin 200 201 write(6,1230) 202 write(6,1241) lattice_unita(1,1), 203 > lattice_unita(2,1), 204 > lattice_unita(3,1) 205 write(6,1242) lattice_unita(1,2), 206 > lattice_unita(2,2), 207 > lattice_unita(3,2) 208 write(6,1243) lattice_unita(1,3), 209 > lattice_unita(2,3), 210 > lattice_unita(3,3) 211 write(6,1244) lattice_unitg(1,1), 212 > lattice_unitg(2,1), 213 > lattice_unitg(3,1) 214 write(6,1245) lattice_unitg(1,2), 215 > lattice_unitg(2,2), 216 > lattice_unitg(3,2) 217 write(6,1246) lattice_unitg(1,3), 218 > lattice_unitg(2,3), 219 > lattice_unitg(3,3) 220 write(6,1231) lattice_omega() 221 write(6,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 222 > pack_nwave_all(0),pack_nwave(0) 223 write(6,1251) lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 224 > pack_nwave_all(1),pack_nwave(1) 225 write(6,*) 226 write(6,*) 227 call util_flush(6) 228 end if 229 230 231 !**** translate system if origin is not zero **** 232 if ((origin(1).ne. 0.0d0).or. 233 > (origin(2).ne. 0.0d0).or . 234 > (origin(3).ne. 0.0d0) ) then 235 236 if (taskid.eq.MASTER) then 237 write(*,*) "...translating origin..." 238 write(*,*) 239 write(*,*) 240 end if 241 242 origin(1) = -origin(1) !* translate system by -origin * 243 origin(2) = -origin(2) 244 origin(3) = -origin(3) 245 246 call ion_translate(origin) 247 call psi_translate(origin, 248 > npack1,(ne(1)+ne(2)), 249 > dcpl_mb(psi2(1))) 250 251 call phafac() !** recompute phase factors ** 252 253 end if 254 255 256 call dplot_gen_psi_dn(ispin,ne, 257 > npack1,nfft3d,nemax, 258 > dcpl_mb(psi2(1)), 259 > dbl_mb(dn(1)), 260 > dcpl_mb(psir(1)), 261 > fractional,dbl_mb(occ2(1))) 262 263c used to localize psi....input needs to be generated for this capabilitiy 264c call psi_dmatrix_localize(ispin,ne,ne,n2ft3d, 265c > dcpl_mb(psir(1)),dcpl_mb(psir2(1))) 266 267 call dplot_loop(rtdb, 268 > ispin,ne, 269 > npack1,nfft3d,nemax, 270 > dcpl_mb(psi2(1)), 271 > dbl_mb(dn(1)), 272 > dcpl_mb(psir(1)), 273 > .false.,tag) 274 275 value = value.and.btdb_get(rtdb,'pspw_dplot:count', 276 > mt_int,1,count) 277 278 279* **** deallocate heap memory **** 280 call strfac_end() 281 if (control_version().eq.3) call coulomb_end() 282 if (control_version().eq.4) call coulomb2_end() 283 call ke_end() 284 call ion_write(rtdb) !*** can also use call ion_destroy()???? 285 !call ion_destroy(rtdb) 286 call ion_end() 287 call mask_end() 288 call Pack_end() 289 call G_end() 290 291 value = BA_free_heap(psir2(2)) 292 value = value.and.BA_free_heap(psir(2)) 293 value = value.and.BA_free_heap(dn(2)) 294 value = value.and.BA_free_heap(psi2(2)) 295 if (fractional) then 296 value = BA_free_heap(occ2(2)) 297 end if 298 299 call D3dB_end(1) 300 if (control_version().eq.4) call D3dB_end(2) 301 call Dne_end() 302 303 if (.not. value) call errquit('pspw_dplot:freeing heap memory',0, 304 & MA_ERR) 305 306 307 if (taskid.eq.MASTER) then 308 call current_second(cpu2) 309 write(6,*) 310 write(6,*) '-----------------' 311 write(6,*) 'cputime in seconds' 312 write(6,*) 'total : ',(cpu2-cpu1) 313 write(6,*) 314 call nwpw_message(4) 315 end if 316 call Parallel_Finalize() 317 318 pspw_dplot = value 319 320 return 321 322*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 323 1000 FORMAT(10X,'****************************************************') 324 1010 FORMAT(10X,'* *') 325 1020 FORMAT(10X,'* pspw DPLOT *') 326 1030 FORMAT(10x,'* [ Generates density and orbital grids ] *') 327 1035 FORMAT(10x,'* [ NorthWest Chemistry implementation ] *') 328 1040 FORMAT(10X,'* version #1.00 08/22/01 *') 329 1100 FORMAT(//) 330 1110 FORMAT(10X,'============ PSPW DPLOT input data =================') 331 1111 FORMAT(/' number of processors used:',I3) 332 1112 FORMAT( ' parallel mapping : 1d slab') 333 1113 FORMAT( ' parallel mapping :2d hilbert') 334 1114 FORMAT( ' parallel mapping : balanced') 335 1115 FORMAT(/' options:') 336 1116 FORMAT( ' parallel mapping : not balanced') 337 1118 FORMAT( ' parallel mapping : 2d hcurve') 338 1121 FORMAT(5X,' boundary conditions = ',A,'(version', I1,')') 339 1130 FORMAT(5X,' electron spin = ',A) 340 1220 FORMAT(/' number of electrons: spin up=',I3,' spin down=',I3,A) 341 1224 FORMAT(/' ncell = ',3I2) 342 1225 FORMAT(/' position tolerance = ',E12.6) 343 1226 FORMAT(/5x,' origin=<',3f8.3,' >') 344 345 1230 FORMAT(/' supercell:') 346 1231 FORMAT(5x,' volume : ',F10.1) 347 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 348 1242 FORMAT(5x,' a2=<',3f8.3,' >') 349 1243 FORMAT(5x,' a3=<',3f8.3,' >') 350 1244 FORMAT(5x,' b1=<',3f8.3,' >') 351 1245 FORMAT(5x,' b2=<',3f8.3,' >') 352 1246 FORMAT(5x,' b3=<',3f8.3,' >') 353 354 1250 FORMAT(5X,' density cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 355 & '( ',I8,' waves ',I8,' per task)') 356 1251 FORMAT(5X,' wavefnc cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 357 & '( ',I8,' waves ',I8,' per task)') 358 359 end 360 361* *********************************************** 362* * * 363* * dplot_gen_psi_dn * 364* * * 365* *********************************************** 366 367 subroutine dplot_gen_psi_dn(ispin,ne, 368 > npack1,nfft3d,nemax, 369 > psi, 370 > dn, 371 > psi_r, 372 > fractional,occ) 373 implicit none 374 integer ispin,ne(2) 375 integer npack1,nfft3d,nemax 376 complex*16 psi(npack1,nemax) 377 real*8 dn(2*nfft3d,2) 378 real*8 psi_r(2*nfft3d,nemax) 379 logical fractional 380 real*8 occ(*) 381 382 383* **** local variables **** 384 integer taskid 385 integer MASTER 386 parameter (MASTER=0) 387 388 integer n2ft3d 389 integer i,ms,n 390 integer n1(2),n2(2) 391 integer nx,ny,nz 392 real*8 scal1,scal2 393 394 395* **** external functions **** 396 integer control_version 397 real*8 lattice_omega 398 external control_version 399 external lattice_omega 400 401 402 call Parallel_taskid(taskid) 403 n2ft3d = 2*nfft3d 404 405 n1(1) = 1 406 n2(1) = ne(1) 407 n1(2) = ne(1) + 1 408 n2(2) = ne(1) + ne(2) 409 410 call D3dB_nx(1,nx) 411 call D3dB_ny(1,ny) 412 call D3dB_nz(1,nz) 413 scal1 = 1.0d0/dble(nx*ny*nz) 414 scal2 = 1.0d0/lattice_omega() 415 416* ******************* 417* **** get psi_r **** 418* ******************* 419 do n=n1(1),n2(ispin) 420 call Pack_c_Copy(1,psi(1,n),psi_r(1,n)) 421 call Pack_c_unpack(1,psi_r(1,n)) 422 call D3dB_cr_fft3b(1,psi_r(1,n)) 423 call D3dB_r_Zero_Ends(1,psi_r(1,n)) 424c call D3dB_r_SMul(1,dsqrt(scal2),psi_r(1,n),psi_r(1,n)) 425 call D3dB_r_SMul1(1,dsqrt(scal2),psi_r(1,n)) 426 end do 427 428 429* ********************* 430* **** generate dn **** 431* ********************* 432 call dcopy(ispin*n2ft3d,0.0d0,0,dn,1) 433 if (fractional) then 434 do ms=1,ispin 435 do n=n1(ms),n2(ms) 436 do i=1,n2ft3d 437c dn(i,ms) = dn(i,ms) + scal2*(psi_r(i,n)**2) 438 dn(i,ms) = dn(i,ms) + (psi_r(i,n)**2)*occ(n) 439 end do 440 end do 441 call D3dB_r_Zero_Ends(1,dn(1,ms)) 442 end do 443 else 444 do ms=1,ispin 445 do n=n1(ms),n2(ms) 446 do i=1,n2ft3d 447c dn(i,ms) = dn(i,ms) + scal2*(psi_r(i,n)**2) 448 dn(i,ms) = dn(i,ms) + (psi_r(i,n)**2) 449 end do 450 end do 451 call D3dB_r_Zero_Ends(1,dn(1,ms)) 452 end do 453 end if 454 455 return 456 end 457 458* *********************************************** 459* * * 460* * dplot_loop * 461* * * 462* *********************************************** 463 464 subroutine dplot_loop(rtdb, 465 > ispin,ne, 466 > npack1,nfft3d,nemax, 467 > psi, 468 > dn, 469 > psi_r, 470 > add_tag,tag) 471 implicit none 472 integer rtdb 473 integer ispin,ne(2) 474 integer npack1,nfft3d,nemax 475 complex*16 psi(npack1,nemax) 476 real*8 dn(2*nfft3d,2) 477 real*8 psi_r(2*nfft3d,nemax) 478 logical add_tag 479 character*8 tag 480 481#include "bafdecls.fh" 482#include "btdb.fh" 483#include "errquit.fh" 484 485* **** local variables **** 486 logical value,grid3d,grid1d 487 integer taskid 488 integer MASTER 489 parameter (MASTER=0) 490 491 integer n2ft3d 492 integer i,i1,i2,i3,count,number,ia,number1,number2,ms 493 integer n1(2),n2(2),nn(2) 494 integer nx,ny,nz,indx,indxr 495 integer rho(2),rhog(2) 496 real*8 scal1,scal2,x(3),r,c 497 498 character*72 cube_comment 499 character*50 name1,name2 500 character*50 filename 501 integer name1_len,name2_len,ind,ind2 502 503 integer n_truncate,idx_truncate(2),rgrid(2),trunc(2) 504 505* **** external functions **** 506 integer control_version 507 real*8 lattice_omega 508 external control_version 509 external lattice_omega 510 511 512 grid3d = .false. 513 if (btdb_get(rtdb,'pspw_dplot:3d_grid:nx',mt_int,1,i)) 514 > grid3d = .true. 515 516 grid1d = .false. 517 if (btdb_get(rtdb,'pspw_dplot:1d_grid:nx',mt_int,1,i)) 518 > grid1d = .true. 519 520 call Parallel_taskid(taskid) 521 n2ft3d = 2*nfft3d 522 value = BA_push_get(mt_dbl,(n2ft3d),'rho',rho(2),rho(1)) 523 if (.not. value) call errquit('out of stack memory',0, MA_ERR) 524 525 n1(1) = 1 526 n2(1) = ne(1) 527 n1(2) = ne(1) + 1 528 n2(2) = ne(1) + ne(2) 529 530 call D3dB_nx(1,nx) 531 call D3dB_ny(1,ny) 532 call D3dB_nz(1,nz) 533 scal1 = 1.0d0/dble(nx*ny*nz) 534 scal2 = 1.0d0/lattice_omega() 535 536* **** setup atom_truncate **** 537 n_truncate = 0 538 if (btdb_get(rtdb,'pspw_dplot:atom_truncate_size', 539 > mt_int,1,n_truncate)) then 540 if (n_truncate.gt.0) then 541 if (BA_push_get(mt_dbl,n2ft3d,'trunc',trunc(2),trunc(1))) then 542 if (BA_push_get(mt_int,n_truncate,'idx_truncate', 543 > idx_truncate(2),idx_truncate(1))) then 544 if (.not.btdb_get(rtdb,'pspw_dplot:atom_truncate', 545 > mt_int,n_truncate,int_mb(idx_truncate(1)))) then 546 n_truncate = 0 547 else 548 if (.not.BA_push_get(mt_dbl,3*n2ft3d, 549 > 'rgrd',rgrid(2),rgrid(1))) 550 > call errquit('pspw_dplot_loop:push stack',0,MA_ERR) 551 call lattice_r_grid(dbl_mb(rgrid(1))) 552 call Truncating_Function_init(rtdb) 553 call Truncating_Function_index(n_truncate, 554 > int_mb(idx_truncate(1)),dbl_mb(rgrid(1)), 555 > dbl_mb(trunc(1))) 556 call Truncating_Function_end() 557 if (.not.BA_pop_stack(rgrid(2))) 558 > call errquit('pspw_dplot_loop:pop stack',0,MA_ERR) 559 end if 560 if (.not.BA_pop_stack(idx_truncate(2))) 561 > call errquit('error popping idx_truncate',0,MA_ERR) 562 end if 563 end if 564 end if 565 end if 566 567 568* ******************************************** 569* **** loop over orbital and density list **** 570* ******************************************** 571 572 value = value.and.btdb_get(rtdb,'pspw_dplot:count', 573 > mt_int,1,count) 574 do i=1,count 575 576* **** define name - not very elegent and could break if **** 577* **** count becomes very large **** 578 ia = ICHAR('a') 579 name1 = 'pspw_dplot:filename'//CHAR(i-1+ia) 580 name2 = 'pspw_dplot:number'//CHAR(i-1+ia) 581 name1_len = index(name1,' ') - 1 582 name2_len = index(name2,' ') - 1 583 584 value = btdb_cget(rtdb,name1(1:name1_len),1,filename) 585 ind = index(filename,' ') - 1 586 if (add_tag) then 587 ind2 = index(tag,' ') - 1 588 filename = filename(1:ind)//tag(1:ind2)//'.cube' 589 ind = index(filename,' ') - 1 590 end if 591 592 if(.not.btdb_get(rtdb,name2(1:name2_len),mt_int,1,number)) then 593 number = 99999999 594 name1 = 'pspw_dplot:number1'//CHAR(i-1+ia) 595 name2 = 'pspw_dplot:number2'//CHAR(i-1+ia) 596 name1_len = index(name1,' ') - 1 597 name2_len = index(name2,' ') - 1 598 value=value.and.btdb_get(rtdb,name1(1:name1_len),mt_int,1, 599 > number1) 600 value=value.and.btdb_get(rtdb,name2(1:name2_len),mt_int,1, 601 > number2) 602 end if 603 if (.not.value) 604 > call errquit( 605 > 'pspw_dplot: btdb_get failed for orbital', 0, RTDB_ERR) 606 607* **** outputing density **** 608 if (number.lt.0) then 609 number = -number 610 611 goto (710,720,730,740,750,760,770,780,785,786 ) number 612 call errquit( 613 > 'dplot_loop: unimplemented directive', number, INPUT_ERR) 614 615* ************* 616* *** total *** 617* ************* 618 710 call D3dB_rr_Sum(1,dn(1,1),dn(1,ispin),dbl_mb(rho(1))) 619 if (.not.add_tag) then 620 if (taskid.eq.MASTER) then 621 write(*,*) ' writing total density', 622 > ' to filename: ',filename(1:ind) 623 end if 624 end if 625 cube_comment = "SCF Total Density" 626 goto 790 627 628* ****************** 629* *** difference *** 630* ****************** 631 720 call D3dB_rr_Sub(1,dn(1,1),dn(1,ispin),dbl_mb(rho(1))) 632 if (.not.add_tag) then 633 if (taskid.eq.MASTER) then 634 write(*,*) ' writing difference density', 635 > ' to filename: ',filename(1:ind) 636 end if 637 end if 638 cube_comment = "SCF Spin Density" 639 goto 790 640 641* ************* 642* *** alpha *** 643* ************* 644 730 call D3dB_r_Copy(1,dn(1,1),dbl_mb(rho(1))) 645 if (.not.add_tag) then 646 if (taskid.eq.MASTER) then 647 write(*,*) ' writing alpha density', 648 > ' to filename: ',filename(1:ind) 649 end if 650 end if 651 cube_comment = "SCF Alpha Density" 652 goto 790 653 654* ************ 655* *** beta *** 656* ************ 657 740 call D3dB_r_Copy(1,dn(1,ispin),dbl_mb(rho(1))) 658 if (.not.add_tag) then 659 if (taskid.eq.MASTER) then 660 write(*,*) ' writing beta density', 661 > ' to filename: ',filename(1:ind) 662 end if 663 end if 664 cube_comment = "SCF Beta Density" 665 goto 790 666 667* ***************** 668* *** laplacian *** 669* ***************** 670 750 call D3dB_rr_Sum(1,dn(1,1),dn(1,ispin),dbl_mb(rho(1))) 671c call D3dB_r_SMul(1,scal1,dbl_mb(rho(1)),dbl_mb(rho(1))) 672 call D3dB_r_SMul1(1,scal1,dbl_mb(rho(1))) 673 call D3dB_rc_fft3f(1,dbl_mb(rho(1))) 674 call Pack_c_pack(1,dbl_mb(rho(1))) 675 nn(1) = 1 676 nn(2) = 0 677 call ke(1,nn,dbl_mb(rho(1)),dbl_mb(rho(1))) 678 call Pack_c_unpack(1,dbl_mb(rho(1))) 679 call D3dB_cr_fft3b(1,dbl_mb(rho(1))) 680 if (.not.add_tag) then 681 if (taskid.eq.MASTER) then 682 write(*,*) ' writing laplacian density', 683 > ' to filename: ',filename(1:ind) 684 end if 685 end if 686 cube_comment = "SCF Laplacian Density" 687 goto 790 688 689* ******************************* 690* *** Electrostatic Potential *** 691* ******************************* 692 760 call D3dB_rr_Sum(1,dn(1,1),dn(1,ispin),dbl_mb(rho(1))) 693 if (control_version().eq.3) then 694 call D3dB_rc_fft3f(1,dbl_mb(rho(1))) 695 call Pack_c_pack(0,dbl_mb(rho(1))) 696c call Pack_c_SMul(0,scal1,dbl_mb(rho(1)),dbl_mb(rho(1))) 697 call Pack_c_SMul1(0,scal1,dbl_mb(rho(1))) 698 call pspw_add_core_dng(1.0d0,dbl_mb(rho(1))) 699c call Pack_c_SMul(0,scal2,dbl_mb(rho(1)),dbl_mb(rho(1))) 700 call Pack_c_SMul1(0,scal2,dbl_mb(rho(1))) 701 call coulomb_v(dbl_mb(rho(1)),dbl_mb(rho(1))) 702c call Pack_c_SMul(0,(1.0d0/scal2), 703c > dbl_mb(rho(1)),dbl_mb(rho(1))) 704 call Pack_c_unpack(0,dbl_mb(rho(1))) 705 call D3dB_cr_fft3b(1,dbl_mb(rho(1))) 706 else 707 call coulomb2_v(dbl_mb(rho(1)),dbl_mb(rho(1))) 708 call pspw_add_core_pot(1.0d0,dbl_mb(rho(1))) 709 end if 710 if (.not.add_tag) then 711 if (taskid.eq.MASTER) then 712 write(*,*) ' writing electrostatic potential', 713 > ' to filename: ',filename(1:ind) 714 end if 715 end if 716 cube_comment = "SCF Electrostatic Potential" 717 goto 790 718 719* ******************** 720* *** ELF *** 721* ******************** 722 770 call generate_ELF(npack1,ne(1), 723 > psi,dn,dbl_mb(rho(1))) 724 if (.not.add_tag) then 725 if (taskid.eq.MASTER) then 726 write(*,*) ' writing restricted/up spin ELF', 727 > ' to filename: ',filename(1:ind) 728 end if 729 end if 730 cube_comment = "SCF ELF" 731 goto 790 732 733 780 call generate_ELF(npack1,ne(2), 734 > psi(1,ne(1)+1),dn(1,2),dbl_mb(rho(1))) 735 if (.not.add_tag) then 736 if (taskid.eq.MASTER) then 737 write(*,*) ' writing down spin ELF', 738 > ' to filename: ',filename(1:ind) 739 end if 740 end if 741 cube_comment = "SCF ELF" 742 goto 790 743 744* *************************** 745* *** density matrix *** 746* *************************** 747 785 name1 = 'pspw_dplot:dmatrix_ms'//CHAR(i-1+ia) 748 name2 = 'pspw_dplot:dmatrix_xyz'//CHAR(i-1+ia) 749 name1_len = index(name1,' ') - 1 750 name2_len = index(name2,' ') - 1 751 if (.not.btdb_get(rtdb,name1(1:name1_len),mt_int,1,ms)) 752 > ms = 1 753 if (.not.btdb_get(rtdb,name2(1:name2_len),mt_dbl,3,x)) then 754 x(1) = 0.0d0 755 x(2) = 0.0d0 756 x(3) = 0.0d0 757 end if 758 if (.not.add_tag) then 759 if (taskid.eq.MASTER) then 760 write(*,*) ' writing density_matrix(ms,:,(x,y,z)', 761 > ' to filename: ',filename(1:ind) 762 write(*,*) "ms =",ms 763 write(*,'(A,3F10.3)') "x,y,z=",x 764 end if 765 end if 766 call generate_dmatrix_column(ms,x, 767 > ispin,ne,n2ft3d, 768 > psi_r,dbl_mb(rho(1))) 769 cube_comment = 'SCF density matrix' 770 goto 790 771 772* *************************** 773* *** structure factor *** 774* *************************** 775 786 if(.not.BA_push_get(mt_dcpl,(nfft3d),'rhog',rhog(2),rhog(1))) 776 > call errquit('out of stack memory',0, MA_ERR) 777 778 call D3dB_rr_Sum(1,dn(1,1),dn(1,ispin),dcpl_mb(rhog(1))) 779 call D3dB_rc_fft3f(1,dcpl_mb(rhog(1))) 780 do i3=1,nz 781 do i2=1,ny 782 do i1=1,(nx/2+1) 783 indx = (i1-1) + (nx/2+1)*(i2-1) + (nx/2+1)*ny*(i3-1) 784 r = dble(dcpl_mb(rhog(1)+indx)) 785 c = dimag(dcpl_mb(rhog(1)+indx)) 786 r = dsqrt(r*r + c*c) 787 indxr = (i1-1) + (nx+2)*(i2-1) + (nx+2)*ny*(i3-1) 788 dbl_mb(rho(1)+indxr) = r 789 if (i>1) then 790 indxr = (nx-i1) + (nx+2)*(i2-1) + (nx+2)*ny*(i3-1) 791 dbl_mb(rho(1)+indxr) = r 792 end if 793 end do 794 end do 795 end do 796 cube_comment = 'Structure Factor' 797 798 if(.not.BA_pop_stack(rhog(2))) 799 > call errquit('popping stack memory',0, MA_ERR) 800 801 goto 790 802 803 804 790 continue 805 number = -number 806 807* **** outputing wavefunction or wavefunction2 **** 808 else 809 !*** used to get spin down homo or lumo **** 810 if (number.eq.123456789) number = ne(1)+1 811 812 !*** overlap functions **** 813 if (number.eq.99999999) then 814 number = -99 815 if ((number1.gt.nemax).or.(number1.lt.1)) then 816 if (taskid.eq.MASTER) then 817 write(*,*) ' Bad orbital number1 ', number1, 818 > ', changing to orbital number1 1.' 819 end if 820 number1 = 1 821 end if 822 if ((number2.gt.nemax).or.(number2.lt.1)) then 823 if (taskid.eq.MASTER) then 824 write(*,*) ' Bad orbital number2 ', number2, 825 > ', changing to orbital number2 1.' 826 end if 827 number2 = 1 828 end if 829 call D3dB_rr_Mul(1,psi_r(1,number1),psi_r(1,number2), 830 > dbl_mb(rho(1))) 831 if (.not.add_tag) then 832 if (taskid.eq.MASTER) then 833 write(*,*) ' writing orbital2 ',number1, number2, 834 > ' to filename: ',filename(1:ind) 835 end if 836 end if 837 cube_comment = 'SCF Molecular Orbitals squared' 838 839 !*** regular molecular orbital **** 840 else 841 if ((number.gt.nemax).or.(number.lt.1)) then 842 if (taskid.eq.MASTER) then 843 write(*,*) ' Bad orbital number ', number, 844 > ', changing to orbital number 1.' 845 end if 846 number = 1 847 end if 848 849 call D3dB_r_Copy(1,psi_r(1,number),dbl_mb(rho(1))) 850 if (.not.add_tag) then 851 if (taskid.eq.MASTER) then 852 write(*,*) ' writing orbital ',number, 853 > ' to filename: ',filename(1:ind) 854 end if 855 end if 856 cube_comment = 'SCF Molecular Orbitals' 857 end if 858 end if 859 860* *** atom_truncate *** 861 if (n_truncate.gt.0) then 862 call D3dB_rr_Mul2(1,dbl_mb(trunc(1)),dbl_mb(rho(1))) 863 end if 864 865 if (grid3d) then 866 call dplot_gcube_write3d(rtdb,filename, 867 > number,cube_comment,dbl_mb(rho(1))) 868 else if (grid1d) then 869 call dplot_gcube_write1d(rtdb,filename, 870 > number,cube_comment,dbl_mb(rho(1))) 871 else 872 call dplot_gcube_write(rtdb,filename, 873 > number,cube_comment,dbl_mb(rho(1))) 874 endif 875 876 end do 877 878* **** dealocate MA local variables **** 879 if (n_truncate.gt.0) then 880 if (.not.BA_pop_stack(trunc(2))) 881 > call errquit('error popping trunc',0,MA_ERR) 882 end if 883 884 value = BA_pop_stack(rho(2)) 885 if (.not. value) call errquit('popping of stack memory',0, MA_ERR) 886 887 888 return 889 end 890 891 892 893* *********************************************** 894* * * 895* * dplot_gcube_write * 896* * * 897* *********************************************** 898 899 subroutine dplot_gcube_write(rtdb,filename, 900 > number,cube_comment,rho) 901 implicit none 902 integer rtdb 903 character*50 filename 904 integer number 905 character*72 cube_comment 906 real*8 rho(*) 907 908 909#include "bafdecls.fh" 910#include "btdb.fh" 911#include "geom.fh" 912#include "errquit.fh" 913 914* **** local variables **** 915 logical value 916 integer unit_id 917 parameter (unit_id = 72) 918 919 integer taskid 920 integer MASTER 921 parameter (MASTER=0) 922 923 924* **** geometry variables **** 925c integer geom1 926c character*16 t 927 integer nion,nion2,ncell(3),nl(3),nh(3),n1,n2,n3 928 real*8 q,rxyz(3) 929 real*8 position_tolerance 930 931* **** lattice variables **** 932 integer np1,np2,np3 933 real*8 ua(3,3),r0(3) 934 935 integer l,orb_flag,i 936 character*255 full_filename 937 integer nfft3d,tmp1(2),tmp2(2) 938 939* **** external functions **** 940 integer ion_nion,control_mapping 941 real*8 lattice_unita,ion_rion,ion_q 942 external ion_nion,control_mapping 943 external lattice_unita,ion_rion,ion_q 944 945 call Parallel_taskid(taskid) 946 947 948* **** OPEN cube FILE **** 949 if (taskid.eq.MASTER) then 950 call util_file_name_noprefix(filename,.false., 951 > .false., 952 > full_filename) 953 l = index(full_filename,' ') -1 954 OPEN(unit_id,file=full_filename(1:l),form='formatted') 955 WRITE(unit_id,*) 'molecule' 956 WRITE(unit_id,*) cube_comment 957 end if 958 959* **** open geom **** 960 nion = ion_nion() 961c value = geom_create(geom1,'geometry') 962c value = value.and.geom_rtdb_load(rtdb,geom1,'geometry') 963c value = value.and.geom_ncent(geom1,nion) 964c if (.not. value) call errquit('opening geometry',0) 965 966* **** write lattice **** 967 call D3dB_nx(1,np1) 968 call D3dB_ny(1,np2) 969 call D3dB_nz(1,np3) 970 do i=1,3 971 ua(i,1) = lattice_unita(i,1)/np1 972 ua(i,2) = lattice_unita(i,2)/np2 973 ua(i,3) = lattice_unita(i,3)/np3 974 975 r0(i) = -( lattice_unita(i,1) 976 > + lattice_unita(i,2) 977 > + lattice_unita(i,3) )/2.0d0 978 end do 979 980* **** get position_tolerance and nion2 = nion+special_positions **** 981 if (.not. 982 > btdb_get(rtdb,'pspw_dplot:position_tolerance', 983 > mt_dbl,1,position_tolerance)) 984 > position_tolerance = 0.0d0 985 986* **** get ncell **** 987 if (.not.btdb_get(rtdb,'pspw_dplot:ncell',mt_int,3,ncell)) then 988 ncell(1) = 0 989 ncell(2) = 0 990 ncell(3) = 0 991 end if 992 if (ncell(1).eq.0) then 993 nl(1) = 0 994 nh(1) = 0 995 else if (ncell(1).gt.0) then 996 nl(1) = -ncell(1) 997 nh(1) = ncell(1) 998 else 999 nl(1) = 0 1000 nh(1) = -ncell(1) 1001 end if 1002 if (ncell(2).eq.0) then 1003 nl(2) = 0 1004 nh(2) = 0 1005 else if (ncell(2).gt.0) then 1006 nl(2) = -ncell(2) 1007 nh(2) = ncell(2) 1008 else 1009 nl(2) = 0 1010 nh(2) = -ncell(2) 1011 end if 1012 if (ncell(3).eq.0) then 1013 nl(3) = 0 1014 nh(3) = 0 1015 else if (ncell(3).gt.0) then 1016 nl(3) = -ncell(3) 1017 nh(3) = ncell(3) 1018 else 1019 nl(3) = 0 1020 nh(3) = -ncell(3) 1021 end if 1022 nion2 = 0 1023 do n3=nl(3),nh(3) 1024 do n2=nl(2),nh(2) 1025 do n1=nl(1),nh(1) 1026 nion2= nion2+nion 1027 end do 1028 end do 1029 end do 1030 1031 1032c nion2 = nion 1033c do i=1,nion 1034cc value = geom_cent_get(geom1,i,t,rxyz,q) 1035cc if (.not. value) call errquit('reading geometry',0) 1036c rxyz(1) = ion_rion(1,i) 1037c rxyz(2) = ion_rion(2,i) 1038c rxyz(3) = ion_rion(3,i) 1039c call special_position_count(position_tolerance, 1040c > rxyz,nion2) 1041c end do 1042 1043 if (taskid.eq.MASTER) then 1044 orb_flag = 1 1045 if (number.gt.0) orb_flag = -1 1046 1047 write(unit_id,'(I5,3F12.6)') orb_flag*nion2,r0 1048 write(unit_id,'(I5,3F12.6)') np1,ua(1,1),ua(2,1),ua(3,1) 1049 write(unit_id,'(I5,3F12.6)') np2,ua(1,2),ua(2,2),ua(3,2) 1050 write(unit_id,'(I5,3F12.6)') np3,ua(1,3),ua(2,3),ua(3,3) 1051 end if 1052 1053 1054 1055 1056* **** write geometry **** 1057 do i=1,nion 1058c value = geom_cent_get(geom1,i,t,rxyz,q) 1059c if (.not. value) call errquit('reading geometry',0) 1060 q = ion_q(i) 1061 do n3=nl(3),nh(3) 1062 do n2=nl(2),nh(2) 1063 do n1=nl(1),nh(1) 1064 rxyz(1) = ion_rion(1,i) + n1*lattice_unita(1,1) 1065 > + n2*lattice_unita(1,2) 1066 > + n3*lattice_unita(1,3) 1067 rxyz(2) = ion_rion(2,i) + n1*lattice_unita(2,1) 1068 > + n2*lattice_unita(2,2) 1069 > + n3*lattice_unita(2,3) 1070 rxyz(3) = ion_rion(3,i) + n1*lattice_unita(3,1) 1071 > + n2*lattice_unita(3,2) 1072 > + n3*lattice_unita(3,3) 1073 if (taskid.eq.MASTER) then 1074 WRITE(unit_id,'(I5,4F12.6)') int(q),q,rxyz 1075 end if 1076 end do 1077 end do 1078 end do 1079 1080c if (taskid.eq.MASTER) then 1081c WRITE(unit_id,'(I5,4F12.6)') int(q),q,rxyz 1082c call special_position_tolerance(position_tolerance, 1083c > unit_id,q,rxyz) 1084c end if 1085 end do 1086 1087* **** write orbital header **** 1088 if (number.gt.0) then 1089 if (taskid.eq.MASTER) write(unit_id,*) 1,number 1090 end if 1091 1092* **** allocate space **** 1093 call D3dB_nfft3d(1,nfft3d) 1094 value = BA_push_get(mt_dcpl,(nfft3d), 'ffttmp1',tmp1(2),tmp1(1)) 1095 value = value.and. 1096 > BA_push_get(mt_dcpl,(nfft3d), 'ffttmp2',tmp2(2),tmp2(1)) 1097 if (.not. value) call errquit('out of stack memory',0, MA_ERR) 1098 1099* **** transpose grid **** 1100 if (control_mapping().eq.1) 1101 > call D3dB_c_transpose_jk(1,rho,dcpl_mb(tmp2(1)),dcpl_mb(tmp1(1))) 1102 1103* **** write grid **** 1104 call D3dB_r_FormatWrite_reverse(1,unit_id,rho,dcpl_mb(tmp1(1))) 1105 1106* **** deallocate space **** 1107 value = BA_pop_stack(tmp2(2)) 1108 value = value.and.BA_pop_stack(tmp1(2)) 1109 if (.not. value) call errquit('popping stack memory',0, MA_ERR) 1110 1111 1112 1113* **** close geom **** 1114c value = geom_destroy(geom1) 1115c if (.not. value) call errquit('closing geometry',0) 1116 1117 1118 1119* **** CLOSE cube FILE **** 1120 if (taskid.eq.MASTER) then 1121 CLOSE(unit_id) 1122 end if 1123 1124 return 1125 end 1126 1127* *********************************************** 1128* * * 1129* * special_position_tolerance * 1130* * * 1131* *********************************************** 1132 1133 subroutine special_position_tolerance(position_tolerance, 1134 > unit,q,r2) 1135 implicit none 1136 real*8 position_tolerance 1137 integer unit 1138 real*8 q 1139 real*8 r2(3) 1140 1141* **** Local variables defined **** 1142 real*8 fa1,fa2,fa3 1143 real*8 a(3,3),b(3,3),volume 1144 integer i,j 1145 real*8 rxyz(3) 1146 1147* **** external functions **** 1148 real*8 lattice_unita 1149 external lattice_unita 1150 1151* ***** Determine the unit lattice vectors and distances ****** 1152 do j=1,3 1153 do i=1,3 1154 a(i,j) = lattice_unita(i,j) 1155 end do 1156 end do 1157 1158 b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3) 1159 b(2,1) = a(3,2)*a(1,3) - a(1,2)*a(3,3) 1160 b(3,1) = a(1,2)*a(2,3) - a(2,2)*a(1,3) 1161 b(1,2) = a(2,3)*a(3,1) - a(3,3)*a(2,1) 1162 b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1) 1163 b(3,2) = a(1,3)*a(2,1) - a(2,3)*a(1,1) 1164 b(1,3) = a(2,1)*a(3,2) - a(3,1)*a(2,2) 1165 b(2,3) = a(3,1)*a(1,2) - a(1,1)*a(3,2) 1166 b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2) 1167 volume = a(1,1)*b(1,1) 1168 > + a(2,1)*b(2,1) 1169 > + a(3,1)*b(3,1) 1170 1171 volume = 1.0d0/volume 1172 call dscal(9,volume,b,1) 1173 1174* *** Break the Ion positions into the a1, a2, and a3 components *** 1175 fa1 = b(1,1) * r2(1) 1176 > + b(2,1) * r2(2) 1177 > + b(3,1) * r2(3) 1178 1179 fa2 = b(1,2) * r2(1) 1180 > + b(2,2) * r2(2) 1181 > + b(3,2) * r2(3) 1182 1183 fa3 = b(1,3) * r2(1) 1184 > + b(2,3) * r2(2) 1185 > + b(3,3) * r2(3) 1186 1187 1188 if ((fa1+position_tolerance) .GT. (0.5d0)) THEN 1189 rxyz(1) = r2(1) - lattice_unita(1,1) 1190 rxyz(2) = r2(2) - lattice_unita(2,1) 1191 rxyz(3) = r2(3) - lattice_unita(3,1) 1192 WRITE(unit,'(I5,4F12.6)') int(q),q,rxyz 1193 end if 1194 1195 if ((fa1-position_tolerance) .LT. (-0.5d0)) THEN 1196 rxyz(1) = r2(1) + lattice_unita(1,1) 1197 rxyz(2) = r2(2) + lattice_unita(2,1) 1198 rxyz(3) = r2(3) + lattice_unita(3,1) 1199 WRITE(unit,'(I5,4F12.6)') int(q),q,rxyz 1200 end if 1201 1202 if ((fa2+position_tolerance) .GT. (0.5d0)) THEN 1203 rxyz(1) = r2(1) - lattice_unita(1,2) 1204 rxyz(2) = r2(2) - lattice_unita(2,2) 1205 rxyz(3) = r2(3) - lattice_unita(3,2) 1206 WRITE(unit,'(I5,4F12.6)') int(q),q,rxyz 1207 end if 1208 1209 if ((fa2-position_tolerance) .LT. (-0.5d0)) THEN 1210 rxyz(1) = r2(1) + lattice_unita(1,2) 1211 rxyz(2) = r2(2) + lattice_unita(2,2) 1212 rxyz(3) = r2(3) + lattice_unita(3,2) 1213 WRITE(unit,'(I5,4F12.6)') int(q),q,rxyz 1214 end if 1215 1216 1217 if ((fa3+position_tolerance) .GT. (0.5d0)) THEN 1218 rxyz(1) = r2(1) - lattice_unita(1,3) 1219 rxyz(2) = r2(2) - lattice_unita(2,3) 1220 rxyz(3) = r2(3) - lattice_unita(3,3) 1221 WRITE(unit,'(I5,4F12.6)') int(q),q,rxyz 1222 end if 1223 1224 if ((fa3-position_tolerance) .LT. (-0.5d0)) THEN 1225 rxyz(1) = r2(1) + lattice_unita(1,3) 1226 rxyz(2) = r2(2) + lattice_unita(2,3) 1227 rxyz(3) = r2(3) + lattice_unita(3,3) 1228 WRITE(unit,'(I5,4F12.6)') int(q),q,rxyz 1229 end if 1230 1231 return 1232 end 1233 1234* *********************************************** 1235* * * 1236* * special_position_count * 1237* * * 1238* *********************************************** 1239 1240 subroutine special_position_count(position_tolerance, 1241 > r2,count) 1242 implicit none 1243 real*8 position_tolerance 1244 real*8 r2(3) 1245 integer count 1246 1247* **** Local variables defined **** 1248 real*8 fa1,fa2,fa3 1249 real*8 a(3,3),b(3,3),volume 1250 integer i,j 1251 1252* **** external functions **** 1253 real*8 lattice_unita 1254 external lattice_unita 1255 1256* ***** Determine the unit lattice vectors and distances ****** 1257 do j=1,3 1258 do i=1,3 1259 a(i,j) = lattice_unita(i,j) 1260 end do 1261 end do 1262 1263 b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3) 1264 b(2,1) = a(3,2)*a(1,3) - a(1,2)*a(3,3) 1265 b(3,1) = a(1,2)*a(2,3) - a(2,2)*a(1,3) 1266 b(1,2) = a(2,3)*a(3,1) - a(3,3)*a(2,1) 1267 b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1) 1268 b(3,2) = a(1,3)*a(2,1) - a(2,3)*a(1,1) 1269 b(1,3) = a(2,1)*a(3,2) - a(3,1)*a(2,2) 1270 b(2,3) = a(3,1)*a(1,2) - a(1,1)*a(3,2) 1271 b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2) 1272 volume = a(1,1)*b(1,1) 1273 > + a(2,1)*b(2,1) 1274 > + a(3,1)*b(3,1) 1275 1276 volume = 1.0d0/volume 1277 call dscal(9,volume,b,1) 1278 1279* *** Break the Ion positions into the a1, a2, and a3 components *** 1280 fa1 = b(1,1) * r2(1) 1281 > + b(2,1) * r2(2) 1282 > + b(3,1) * r2(3) 1283 1284 fa2 = b(1,2) * r2(1) 1285 > + b(2,2) * r2(2) 1286 > + b(3,2) * r2(3) 1287 1288 fa3 = b(1,3) * r2(1) 1289 > + b(2,3) * r2(2) 1290 > + b(3,3) * r2(3) 1291 1292 1293 if ((fa1+position_tolerance) .GT. (0.5d0)) THEN 1294 count = count+1 1295 end if 1296 1297 if ((fa1-position_tolerance) .LT. (-0.5d0)) THEN 1298 count = count + 1 1299 end if 1300 1301 if ((fa2+position_tolerance) .GT. (0.5d0)) THEN 1302 count = count + 1 1303 end if 1304 1305 if ((fa2-position_tolerance) .LT. (-0.5d0)) THEN 1306 count = count + 1 1307 end if 1308 1309 1310 if ((fa3+position_tolerance) .GT. (0.5d0)) THEN 1311 count = count + 1 1312 end if 1313 1314 if ((fa3-position_tolerance) .LT. (-0.5d0)) THEN 1315 count = count + 1 1316 end if 1317 1318 return 1319 end 1320 1321 1322* ******************************** 1323* * * 1324* * pspw_add_core_dng * 1325* * * 1326* ******************************** 1327 subroutine pspw_add_core_dng(rcut,dng) 1328 implicit none 1329#include "errquit.fh" 1330 real*8 rcut 1331 complex*16 dng(*) 1332 1333#include "bafdecls.fh" 1334 1335* *** local variables *** 1336 logical value 1337 integer nfft3d 1338 integer i,ii 1339 integer exi(2),vg(2),G(3) 1340 real*8 w,scal,gg 1341 1342* **** external functions **** 1343 integer G_indx,ion_nion,ion_katm 1344 real*8 ion_zv,lattice_omega 1345 external G_indx,ion_nion,ion_katm 1346 external ion_zv,lattice_omega 1347 1348 1349 call D3dB_nfft3d(1,nfft3d) 1350 value = BA_push_get(mt_dcpl,nfft3d,'exi', exi(2), exi(1)) 1351 value = value.and. 1352 > BA_push_get(mt_dbl,nfft3d,'vg',vg(2),vg(1)) 1353 if (.not. value) call errquit('pspw_add_core_dng:push stack',0, 1354 & MA_ERR) 1355 1356 G(1) = G_indx(1) 1357 G(2) = G_indx(2) 1358 G(3) = G_indx(3) 1359 w = 0.25d0*rcut*rcut 1360 1361 1362 do ii=1,ion_nion() 1363 scal = -ion_zv(ii)/lattice_omega() 1364 1365* *** fourier transform of a gaussian *** 1366 do i=1,nfft3d 1367 gg = ( dbl_mb(G(1)+i-1)*dbl_mb(G(1)+i-1) 1368 > + dbl_mb(G(2)+i-1)*dbl_mb(G(2)+i-1) 1369 > + dbl_mb(G(3)+i-1)*dbl_mb(G(3)+i-1) ) 1370 1371 dbl_mb(vg(1)+i-1) = (scal)*exp(-w*gg) 1372 end do 1373 1374* **** add to dng *** 1375 call strfac(ii,dcpl_mb(exi(1))) 1376 call Pack_c_pack(0,dcpl_mb(exi(1))) 1377 call Pack_t_pack(0,dbl_mb(vg(1))) 1378c call Pack_tc_Mul(0,dbl_mb(vg(1)), 1379c > dcpl_mb(exi(1)), 1380c > dcpl_mb(exi(1))) 1381c call Pack_cc_Sum(0,dng,dcpl_mb(exi(1)),dng) 1382 call Pack_tc_Mul2(0,dbl_mb(vg(1)),dcpl_mb(exi(1))) 1383 call Pack_cc_Sum2(0,dcpl_mb(exi(1)),dng) 1384 1385 end do 1386 1387 value = BA_pop_stack(vg(2)) 1388 value = value.and.BA_pop_stack(exi(2)) 1389 if (.not. value) call errquit('pspw_add_core_dng:pop stack',1, 1390 & MA_ERR) 1391 return 1392 end 1393 1394* ******************************** 1395* * * 1396* * pspw_add_core_pot * 1397* * * 1398* ******************************** 1399 subroutine pspw_add_core_pot(rcut,vh) 1400 implicit none 1401#include "errquit.fh" 1402 real*8 rcut 1403 real*8 vh(*) 1404 1405#include "bafdecls.fh" 1406 1407* *** local variables *** 1408 logical value 1409 integer n2ft3d 1410 integer i,ii 1411 integer rgrid(2) 1412 real*8 c,q,r,xerf,yerf,sqrt_pi 1413 real*8 xii,yii,zii 1414 real*8 xi,yi,zi 1415 1416* **** external functions **** 1417 integer ion_nion,ion_katm 1418 real*8 ion_zv,ion_rion,util_erf 1419 external ion_nion,ion_katm 1420 external ion_zv,ion_rion,util_erf 1421 1422 1423 call D3dB_n2ft3d(1,n2ft3d) 1424 value = BA_push_get(mt_dbl,3*n2ft3d,'rgrid',rgrid(2),rgrid(1)) 1425 if (.not. value) call errquit('pspw_add_core_pot:push stack',0, 1426 & MA_ERR) 1427 1428 sqrt_pi = dsqrt(4.0d0*datan(1.0d0)) 1429 c = 1.0d0/rcut 1430 1431 call lattice_r_grid(dbl_mb(rgrid(1))) 1432 1433 do ii=1,ion_nion() 1434 q = -ion_zv(ii) 1435 xii = ion_rion(1,ii) 1436 yii = ion_rion(2,ii) 1437 zii = ion_rion(3,ii) 1438 1439 do i=1,n2ft3d 1440 xi = dbl_mb(rgrid(1) + 3*(i-1)) 1441 yi = dbl_mb(rgrid(1) + 3*(i-1)+1) 1442 zi = dbl_mb(rgrid(1) + 3*(i-1)+2) 1443 r = dsqrt((xii-xi)**2 + (yii-yi)**2 + (zii-zi)**2) 1444 1445 if (r .gt. 1.0d-15) then 1446 xerf = r*c 1447 yerf = util_erf(xerf) 1448 vh(i) = vh(i) + (q/r)*yerf 1449 else 1450 vh(i) = vh(i) + 2.0d0*q*c/sqrt_pi 1451 end if 1452 end do 1453 1454 end do 1455 1456 value = BA_pop_stack(rgrid(2)) 1457 if (.not. value) call errquit('pspw_add_core_pot:pop stack',1, 1458 & MA_ERR) 1459 return 1460 end 1461 1462 1463* ******************************** 1464* * * 1465* * psi_translate * 1466* * * 1467* ******************************** 1468 subroutine psi_translate(trans,npack1,nemax,psi) 1469 implicit none 1470 real*8 trans(3) 1471 integer npack1,nemax 1472 complex*16 psi(npack1,nemax) 1473 1474#include "bafdecls.fh" 1475#include "errquit.fh" 1476 1477* **** local variables **** 1478 logical value 1479 integer n,nion 1480 integer exi(2),rion(2) 1481 1482* **** external functions **** 1483 integer ion_nion 1484 real*8 lattice_unita 1485 external ion_nion 1486 external lattice_unita 1487 1488 nion = ion_nion() 1489 value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1)) 1490 value = value.and. 1491 > BA_push_get(mt_dbl,3*nion,'rion', rion(2), rion(1)) 1492 if (.not. value) call errquit('psi_translate:push stack',0,0) 1493 call dcopy(3*nion,0.0d0,0,dbl_mb(rion(1)),1) 1494 1495 dbl_mb(rion(1) ) = trans(1) 1496 > - ( lattice_unita(1,1) 1497 > + lattice_unita(1,2) 1498 > + lattice_unita(1,3) )/2.0d0 1499 dbl_mb(rion(1)+1) = trans(2) 1500 > - ( lattice_unita(2,1) 1501 > + lattice_unita(2,2) 1502 > + lattice_unita(2,3) )/2.0d0 1503 dbl_mb(rion(1)+2) = trans(3) 1504 > - ( lattice_unita(3,1) 1505 > + lattice_unita(3,2) 1506 > + lattice_unita(3,3) )/2.0d0 1507 call phafac_rion(dbl_mb(rion(1))) 1508 call strfac_pack(1,1,dcpl_mb(exi(1))) 1509 1510 1511 do n=1,nemax 1512c call Pack_cc_Mul(1,dcpl_mb(exi(1)),psi(1,n),psi(1,n)) 1513 call Pack_cc_Mul2(1,dcpl_mb(exi(1)),psi(1,n)) 1514 !call Pack_cc_conjgMul(1,dcpl_mb(exi(1)),psi(1,n),psi(1,n)) 1515 end do 1516 1517 1518 value = BA_pop_stack(rion(2)) 1519 value = value.and.BA_pop_stack(exi(2)) 1520 if (.not. value) call errquit('psi_translate:pop stack',1,0) 1521 return 1522 end 1523 1524 1525 1526* 1527* *********************************************** 1528* * * 1529* * dplot_gcube_write3d * 1530* * * 1531* *********************************************** 1532 1533 subroutine dplot_gcube_write3d(rtdb,filename, 1534 > number,cube_comment,rho) 1535 implicit none 1536 integer rtdb 1537 character*50 filename 1538 integer number 1539 character*72 cube_comment 1540 real*8 rho(*) 1541 1542 1543#include "bafdecls.fh" 1544#include "btdb.fh" 1545#include "geom.fh" 1546#include "errquit.fh" 1547 1548* **** local variables **** 1549 logical value 1550 integer unit_id 1551 parameter (unit_id = 72) 1552 1553 integer taskid 1554 integer MASTER 1555 parameter (MASTER=0) 1556 1557 1558* **** geometry variables **** 1559c integer geom1 1560c character*16 t 1561 integer nion,nion2 1562 real*8 q,rxyz(3),p(3) 1563 real*8 position_tolerance 1564 1565* **** lattice variables **** 1566 integer np1,np2,np3 1567 real*8 ua(3,3),r0(3) 1568 1569 integer nx,ny,nz 1570 integer nxh,nyh,nzh 1571 real*8 sx,sy,sz,xy,xz,yz,scal1 1572 real*8 dx,dy,dz,w 1573 real*8 sizex(2),sizey(2),sizez(2) 1574 real*8 xaxis(3),yaxis(3),zaxis(3),qq(3) 1575 1576 1577 integer l,orb_flag,i,j,k,ii,jj,kk,indx 1578 character*255 full_filename 1579 integer nfft3d,tmp1(2),tmp2(2),tmp3(2) 1580 integer Gx,Gy,Gz 1581 1582* **** external functions **** 1583 integer ion_nion 1584 real*8 lattice_unita,ion_rion,ion_q 1585 external ion_nion 1586 external lattice_unita,ion_rion,ion_q 1587 1588 integer G_indx 1589 external G_indx 1590 1591 1592 call Parallel_taskid(taskid) 1593 1594 1595* **** OPEN cube FILE **** 1596 if (taskid.eq.MASTER) then 1597 call util_file_name_noprefix(filename,.false., 1598 > .false., 1599 > full_filename) 1600 l = index(full_filename,' ') -1 1601 OPEN(unit_id,file=full_filename(1:l),form='formatted') 1602 WRITE(unit_id,*) 'molecule' 1603 WRITE(unit_id,*) cube_comment 1604 end if 1605 1606* **** open geom **** 1607 nion = ion_nion() 1608c value = geom_create(geom1,'geometry') 1609c value = value.and.geom_rtdb_load(rtdb,geom1,'geometry') 1610c value = value.and.geom_ncent(geom1,nion) 1611c if (.not. value) call errquit('opening geometry',0) 1612 1613* **** read the origin *** 1614 if (.not. 1615 > btdb_get(rtdb,'pspw_dplot:3d_grid:o', 1616 > mt_dbl,3,r0)) then 1617 r0(1) = 0.0 1618 r0(2) = 0.0 1619 r0(3) = 0.0 1620 end if 1621 1622* **** read the x-axis,y-axis,z-axis **** 1623 if (.not. 1624 > btdb_get(rtdb,'pspw_dplot:3d_grid:x', 1625 > mt_dbl,3,xaxis)) then 1626 xaxis(1) = 1.0 1627 xaxis(2) = 0.0 1628 xaxis(3) = 0.0 1629 end if 1630 if (.not. 1631 > btdb_get(rtdb,'pspw_dplot:3d_grid:y', 1632 > mt_dbl,3,yaxis)) then 1633 yaxis(1) = 0.0 1634 yaxis(2) = 1.0 1635 yaxis(3) = 0.0 1636 end if 1637 if (.not. 1638 > btdb_get(rtdb,'pspw_dplot:3d_grid:z', 1639 > mt_dbl,3,zaxis)) then 1640 zaxis(1) = 0.0 1641 zaxis(2) = 0.0 1642 zaxis(3) = 1.0 1643 end if 1644 1645* ** normalization of x ** 1646 sx = dsqrt(xaxis(1)**2+xaxis(2)**2+xaxis(3)**2) 1647 do i=1,3 1648 xaxis(i) = xaxis(i) / sx 1649 end do 1650 1651* ** orthogonalization of y wrt x ** 1652 xy = xaxis(1)*yaxis(1)+xaxis(2)*yaxis(2)+xaxis(3)*yaxis(3) 1653 do i=1,3 1654 yaxis(i) = yaxis(i) - xy*xaxis(i) 1655 end do 1656 1657* ** normalization of y ** 1658 sy = dsqrt(yaxis(1)**2+yaxis(2)**2+yaxis(3)**2) 1659 do i=1,3 1660 yaxis(i) = yaxis(i) / sy 1661 end do 1662 1663* ** orthogonalization of z wrt x ** 1664 xz = xaxis(1)*zaxis(1)+xaxis(2)*zaxis(2)+xaxis(3)*zaxis(3) 1665 do i=1,3 1666 zaxis(i) = zaxis(i) - xz*xaxis(i) 1667 end do 1668 1669* ** orthogonalization of z wrt y ** 1670 yz = yaxis(1)*zaxis(1)+yaxis(2)*zaxis(2)+yaxis(3)*zaxis(3) 1671 do i=1,3 1672 zaxis(i) = zaxis(i) - yz*yaxis(i) 1673 end do 1674 1675* ** normalization of z ** 1676 sz = dsqrt(zaxis(1)**2+zaxis(2)**2+zaxis(3)**2) 1677 do i=1,3 1678 zaxis(i) = zaxis(i) / sz 1679 end do 1680 1681 1682 1683* **** read nx,ny,and nz **** 1684 if (.not. 1685 > btdb_get(rtdb,'pspw_dplot:3d_grid:nx', 1686 > mt_int,1,nx)) 1687 > nx = 32 1688 if (.not. 1689 > btdb_get(rtdb,'pspw_dplot:3d_grid:ny', 1690 > mt_int,1,ny)) 1691 > ny = 32 1692 if (.not. 1693 > btdb_get(rtdb,'pspw_dplot:3d_grid:nz', 1694 > mt_int,1,nz)) 1695 > nz = 32 1696 nxh = nx/2 1697 nyh = ny/2 1698 nzh = nz/2 1699 1700* **** read sizex,sizey,sizez **** 1701 if (.not. 1702 > btdb_get(rtdb,'pspw_dplot:3d_grid:sizex', 1703 > mt_dbl,2,sizex)) then 1704 sizex(1) = 0.0 1705 sizex(2) = 10.0 1706 end if 1707 if (.not. 1708 > btdb_get(rtdb,'pspw_dplot:3d_grid:sizey', 1709 > mt_dbl,2,sizey)) then 1710 sizey(1) = 0.0 1711 sizey(2) = 10.0 1712 end if 1713 if (.not. 1714 > btdb_get(rtdb,'pspw_dplot:3d_grid:sizez', 1715 > mt_dbl,2,sizez)) then 1716 sizez(1) = 0.0 1717 sizez(2) = 10.0 1718 end if 1719 1720 if (taskid.eq.MASTER) then 1721 write(*,*) 1722 write(*,1240) 1723 write(*,1241) r0 1724 write(*,1242) xaxis,sizex,nx 1725 write(*,1243) yaxis,sizey,ny 1726 write(*,1244) zaxis,sizez,nz 1727 write(*,*) 1728 write(*,*) 1729 end if 1730 1240 FORMAT(5x,' 3d-Grid Generation') 1731 1241 FORMAT(5x,' origin=<',3f8.3,' >') 1732 1242 FORMAT(5x,' x-axis=<',3f8.3,' > xmin,xmax=',2F8.3,' nx=',I3) 1733 1243 FORMAT(5x,' y-axis=<',3f8.3,' > ymin,ymax=',2F8.3,' ny=',I3) 1734 1244 FORMAT(5x,' z-axis=<',3f8.3,' > zmin,zmax=',2F8.3,' nz=',I3) 1735 1736 1737 call D3dB_nx(1,np1) 1738 call D3dB_ny(1,np2) 1739 call D3dB_nz(1,np3) 1740 1741 dx = (sizex(2)-sizex(1))/dble(nx) 1742 dy = (sizey(2)-sizey(1))/dble(ny) 1743 dz = (sizez(2)-sizez(1))/dble(nz) 1744 do i=1,3 1745 ua(i,1) = dx*xaxis(i) 1746 ua(i,2) = dy*yaxis(i) 1747 ua(i,3) = dz*zaxis(i) 1748 end do 1749 1750* ** shift origin ** 1751 qq(1) = r0(1) + 0.5d0*(sizex(1)+sizex(2)) 1752 > - 0.5d0*(lattice_unita(1,1) 1753 > +lattice_unita(1,2) 1754 > +lattice_unita(1,3)) 1755 qq(2) = r0(2) + 0.5d0*(sizey(1)+sizey(2)) 1756 > - 0.5d0*(lattice_unita(2,1) 1757 > +lattice_unita(2,2) 1758 > +lattice_unita(2,3)) 1759 qq(3) = r0(3) + 0.5d0*(sizez(1)+sizez(2)) 1760 > - 0.5d0*(lattice_unita(3,1) 1761 > +lattice_unita(3,2) 1762 > +lattice_unita(3,3)) 1763 1764 r0(1) = r0(1) + sizex(1) 1765 r0(2) = r0(2) + sizey(1) 1766 r0(3) = r0(3) + sizez(1) 1767 1768 1769 1770 1771* **** get position_tolerance and nion2 = nion+special_positions **** 1772 if (.not. 1773 > btdb_get(rtdb,'pspw_dplot:position_tolerance', 1774 > mt_dbl,1,position_tolerance)) 1775 > position_tolerance = 0.0d0 1776 1777 nion2 = nion 1778 do i=1,nion 1779c value = geom_cent_get(geom1,i,t,rxyz,q) 1780c if (.not. value) call errquit('reading geometry',0) 1781 rxyz(1) = ion_rion(1,i) 1782 rxyz(2) = ion_rion(2,i) 1783 rxyz(3) = ion_rion(3,i) 1784 call special_position_count(position_tolerance, 1785 > rxyz,nion2) 1786 end do 1787 1788 1789 if (taskid.eq.MASTER) then 1790 orb_flag = 1 1791 if (number.gt.0) orb_flag = -1 1792 1793 write(unit_id,'(I5,3F12.6)') orb_flag*nion2,r0 1794 write(unit_id,'(I5,3F12.6)') nx,ua(1,1),ua(2,1),ua(3,1) 1795 write(unit_id,'(I5,3F12.6)') ny,ua(1,2),ua(2,2),ua(3,2) 1796 write(unit_id,'(I5,3F12.6)') nz,ua(1,3),ua(2,3),ua(3,3) 1797 end if 1798 1799 1800* **** write geometry **** 1801 do i=1,nion 1802c value = geom_cent_get(geom1,i,t,rxyz,q) 1803c if (.not. value) call errquit('reading geometry',0) 1804 rxyz(1) = ion_rion(1,i) 1805 rxyz(2) = ion_rion(2,i) 1806 rxyz(3) = ion_rion(3,i) 1807 q = ion_q(i) 1808 if (taskid.eq.MASTER) then 1809 WRITE(unit_id,'(I5,4F12.6)') int(q),q,rxyz 1810 call special_position_tolerance(position_tolerance, 1811 > unit_id,q,rxyz) 1812 end if 1813 end do 1814 1815* **** write orbital header **** 1816 if (number.gt.0) then 1817 if (taskid.eq.MASTER) write(unit_id,*) 1,number 1818 end if 1819 1820* **** allocate space **** 1821 call D3dB_nfft3d(1,nfft3d) 1822 value = BA_push_get(mt_dcpl,(nfft3d), 'ffttmp1',tmp1(2),tmp1(1)) 1823 value = value.and. 1824 > BA_push_get(mt_dcpl,(nfft3d), 'ffttmp2',tmp2(2),tmp2(1)) 1825 value = value.and. 1826 > BA_push_get(mt_dbl,(nx*ny*nz),'tmp3',tmp3(2),tmp3(1)) 1827 if (.not. value) call errquit('out of stack memory',0, MA_ERR) 1828 1829 1830* **** fft density **** 1831 scal1=1.0d0/dble(np1*np2*np3) 1832 call D3dB_r_SMul(1,scal1,rho,dcpl_mb(tmp2(1))) 1833 call D3dB_r_Zero_Ends(1,dcpl_mb(tmp2(1))) 1834 call D3dB_rc_fft3f(1,dcpl_mb(tmp2(1))) 1835 call Pack_c_pack(0,dcpl_mb(tmp2(1))) 1836 1837 1838* **** density at mesh points **** 1839 Gx= G_indx(1) 1840 Gy= G_indx(2) 1841 Gz= G_indx(3) 1842 do k=1,nz 1843 do j=1,ny 1844 do i=1,nx 1845 ii = i-1-nxh 1846 jj = j-1-nyh 1847 kk = k-1-nzh 1848 p(1) = ii*ua(1,1) + jj*ua(1,2) + kk*ua(1,3) + qq(1) 1849 p(2) = ii*ua(2,1) + jj*ua(2,2) + kk*ua(2,3) + qq(2) 1850 p(3) = ii*ua(3,1) + jj*ua(3,2) + kk*ua(3,3) + qq(3) 1851 do l=1,nfft3d 1852 w = dbl_mb(Gx+l-1)*p(1) 1853 > + dbl_mb(Gy+l-1)*p(2) 1854 > + dbl_mb(Gz+l-1)*p(3) 1855 dcpl_mb(tmp1(1)+l-1)=dcmplx(dcos(w),-dsin(w)) 1856 end do 1857 call Pack_c_pack(0,dcpl_mb(tmp1(1))) 1858 call Pack_cc_dot(0,dcpl_mb(tmp1(1)),dcpl_mb(tmp2(1)),w) 1859 indx = (i-1) + (j-1)*nx + (k-1)*nx*ny 1860 dbl_mb(tmp3(1)+indx) = w 1861 end do 1862 end do 1863 end do 1864 1865* **** write grid **** 1866 if (taskid.eq.MASTER) then 1867 do i=1,nx 1868 do j=1,ny 1869 write(unit_id,'(6e13.5)') 1870 > (dbl_mb(tmp3(1)+(i-1) 1871 > +(j-1)*nx 1872 > +(k-1)*nx*ny),k=1,nz) 1873 end do 1874 end do 1875 end if 1876 1877* **** deallocate space **** 1878 value = BA_pop_stack(tmp3(2)) 1879 value = value.and.BA_pop_stack(tmp2(2)) 1880 value = value.and.BA_pop_stack(tmp1(2)) 1881 if (.not. value) call errquit('popping stack memory',0, MA_ERR) 1882 1883 1884 1885* **** close geom **** 1886c value = geom_destroy(geom1) 1887c if (.not. value) call errquit('closing geometry',0) 1888 1889 1890 1891* **** CLOSE cube FILE **** 1892 if (taskid.eq.MASTER) then 1893 CLOSE(unit_id) 1894 end if 1895 1896 return 1897 end 1898 1899 1900 1901 1902 1903* 1904* *********************************************** 1905* * * 1906* * dplot_gcube_write1d * 1907* * * 1908* *********************************************** 1909 1910 subroutine dplot_gcube_write1d(rtdb,filename, 1911 > number,cube_comment,rho) 1912 implicit none 1913 integer rtdb 1914 character*50 filename 1915 integer number 1916 character*72 cube_comment 1917 real*8 rho(*) 1918 1919 1920#include "bafdecls.fh" 1921#include "btdb.fh" 1922#include "geom.fh" 1923#include "errquit.fh" 1924 1925* **** local variables **** 1926 logical value 1927 integer unit_id 1928 parameter (unit_id = 72) 1929 1930 integer taskid 1931 integer MASTER 1932 parameter (MASTER=0) 1933 1934 1935* **** geometry variables **** 1936c integer geom1 1937c character*16 t 1938 integer nion,nion2 1939 real*8 q,rxyz(3),p(3) 1940 real*8 position_tolerance 1941 1942* **** lattice variables **** 1943 integer np1,np2,np3 1944 real*8 ua(3,3),r0(3) 1945 1946 integer nx,ny,nz 1947 integer nxh,nyh,nzh 1948 real*8 sx,sy,sz,xy,xz,yz,scal1 1949 real*8 dx,dy,dz,w 1950 real*8 sizex(2),sizey(2),sizez(2) 1951 real*8 xaxis(3),yaxis(3),zaxis(3),qq(3) 1952 1953 1954 integer l,orb_flag,i,j,k,ii,jj,kk,indx 1955 character*255 full_filename 1956 integer nfft3d,tmp1(2),tmp2(2),tmp3(2) 1957 integer Gx,Gy,Gz 1958 1959* **** external functions **** 1960 integer ion_nion 1961 real*8 lattice_unita,ion_rion,ion_q 1962 external ion_nion 1963 external lattice_unita,ion_rion,ion_q 1964 1965 integer G_indx 1966 external G_indx 1967 1968 1969 call Parallel_taskid(taskid) 1970 1971 1972* **** OPEN cube FILE **** 1973 if (taskid.eq.MASTER) then 1974 call util_file_name_noprefix(filename,.false., 1975 > .false., 1976 > full_filename) 1977 l = index(full_filename,' ') -1 1978 OPEN(unit_id,file=full_filename(1:l),form='formatted') 1979c WRITE(unit_id,*) 'molecule' 1980c WRITE(unit_id,*) cube_comment 1981 end if 1982 1983* **** open geom **** 1984 nion = ion_nion() 1985c value = geom_create(geom1,'geometry') 1986c value = value.and.geom_rtdb_load(rtdb,geom1,'geometry') 1987c value = value.and.geom_ncent(geom1,nion) 1988c if (.not. value) call errquit('opening geometry',0) 1989 1990* **** read the origin *** 1991 if (.not. 1992 > btdb_get(rtdb,'pspw_dplot:1d_grid:o', 1993 > mt_dbl,3,r0)) then 1994 r0(1) = 0.0 1995 r0(2) = 0.0 1996 r0(3) = 0.0 1997 end if 1998 1999* **** read the x point **** 2000 if (.not. 2001 > btdb_get(rtdb,'pspw_dplot:1d_grid:x', 2002 > mt_dbl,3,xaxis)) then 2003 xaxis(1) = 1.0 2004 xaxis(2) = 0.0 2005 xaxis(3) = 0.0 2006 end if 2007 2008* **** read nx *** 2009 if (.not. 2010 > btdb_get(rtdb,'pspw_dplot:1d_grid:nx', 2011 > mt_int,1,nx)) 2012 > nx = 32 2013 nxh = nx/2 2014 2015 if (taskid.eq.MASTER) then 2016 write(*,*) 2017 write(*,1240) 2018 write(*,1241) r0 2019 write(*,1242) xaxis,nx 2020 write(*,*) 2021 write(*,*) 2022 end if 2023 1240 FORMAT(5x,' 1d-Grid Generation') 2024 1241 FORMAT(5x,' origin=<',3f8.3,' >') 2025 1242 FORMAT(5x,' x-point=<',3f8.3,' > nx=',I3) 2026 2027* *** define line **** 2028 do i=1,3 2029 xaxis(i) = xaxis(i)-r0(i) 2030 end do 2031 2032* ** normalization of x ** 2033 sx = dsqrt(xaxis(1)**2+xaxis(2)**2+xaxis(3)**2) 2034 do i=1,3 2035 xaxis(i) = xaxis(i) / sx 2036 end do 2037 2038 2039 call D3dB_nx(1,np1) 2040 call D3dB_ny(1,np2) 2041 call D3dB_nz(1,np3) 2042 2043 dx = sx/dble(nx) 2044 do i=1,3 2045 ua(i,1) = dx*xaxis(i) 2046 ua(i,2) = 0.0d0 2047 ua(i,3) = 0.0d0 2048 end do 2049 2050* ** shift origin ** 2051 qq(1) = r0(1) 2052 > - 0.5d0*(lattice_unita(1,1) 2053 > +lattice_unita(1,2) 2054 > +lattice_unita(1,3)) 2055 qq(2) = r0(2) 2056 > - 0.5d0*(lattice_unita(2,1) 2057 > +lattice_unita(2,2) 2058 > +lattice_unita(2,3)) 2059 qq(3) = r0(3) 2060 > - 0.5d0*(lattice_unita(3,1) 2061 > +lattice_unita(3,2) 2062 > +lattice_unita(3,3)) 2063 2064 2065* **** allocate space **** 2066 call D3dB_nfft3d(1,nfft3d) 2067 value = BA_push_get(mt_dcpl,(nfft3d), 'ffttmp1',tmp1(2),tmp1(1)) 2068 value = value.and. 2069 > BA_push_get(mt_dcpl,(nfft3d), 'ffttmp2',tmp2(2),tmp2(1)) 2070 if (.not. value) call errquit('out of stack memory',0, MA_ERR) 2071 2072 2073* **** fft density **** 2074 scal1=1.0d0/dble(np1*np2*np3) 2075 call D3dB_r_SMul(1,scal1,rho,dcpl_mb(tmp2(1))) 2076 call D3dB_r_Zero_Ends(1,dcpl_mb(tmp2(1))) 2077 call D3dB_rc_fft3f(1,dcpl_mb(tmp2(1))) 2078 call Pack_c_pack(0,dcpl_mb(tmp2(1))) 2079 2080 2081* **** density at mesh points **** 2082 Gx= G_indx(1) 2083 Gy= G_indx(2) 2084 Gz= G_indx(3) 2085 do i=1,nx 2086 !ii = i-1-nxh 2087 ii = i-1 2088 p(1) = ii*ua(1,1) + qq(1) 2089 p(2) = ii*ua(2,1) + qq(2) 2090 p(3) = ii*ua(3,1) + qq(3) 2091 do l=1,nfft3d 2092 w = dbl_mb(Gx+l-1)*p(1) 2093 > + dbl_mb(Gy+l-1)*p(2) 2094 > + dbl_mb(Gz+l-1)*p(3) 2095 dcpl_mb(tmp1(1)+l-1)=dcmplx(dcos(w),-dsin(w)) 2096 end do 2097 call Pack_c_pack(0,dcpl_mb(tmp1(1))) 2098 call Pack_cc_dot(0,dcpl_mb(tmp1(1)),dcpl_mb(tmp2(1)),w) 2099 2100 if (taskid.eq.MASTER) 2101 > write(unit_id,'(2e13.5)') (i-1)*dx,w 2102 end do 2103 2104 2105* **** deallocate space **** 2106 value = BA_pop_stack(tmp2(2)) 2107 value = value.and.BA_pop_stack(tmp1(2)) 2108 if (.not. value) call errquit('popping stack memory',0, MA_ERR) 2109 2110 2111 2112* **** CLOSE FILE **** 2113 if (taskid.eq.MASTER) then 2114 CLOSE(unit_id) 2115 end if 2116 2117 return 2118 end 2119 2120 2121