1******************************************************** 2* * 3* Water pseudopotential module * 4* * 5* BLCJ Electron-Water pseudopotential * 6* * 7* Interfaced to nwchem-PSPW code * 8* * 9* -- developed by Eric J. Bylaska on July 19,2001 * 10* * 11********************************************************** 12 13 14* ********************************** 15* * * 16* * Waterpsp_init * 17* * * 18* ********************************** 19 20 subroutine Waterpsp_init(rtdb) 21 implicit none 22#include "errquit.fh" 23 integer rtdb 24 25#include "bafdecls.fh" 26#include "btdb.fh" 27#include "geom.fh" 28 29 30* ***** water common block **** 31#include "waterpsp.fh" 32 33 logical mmexist 34 common / pspwqmm/ mmexist 35 36* **** local variables **** 37 logical value 38 integer i,geom,nion 39 double precision q 40 character*16 t 41 42 integer control_version 43 external control_version 44 45 46* **************************** 47* **** read in water data **** 48* **************************** 49 if (mmexist) then 50 value = geom_create(geom,'qmmmgeometry') 51 value = value.and.geom_rtdb_load(rtdb,geom,'qmmmgeometry') 52 value = value.and.geom_ncent(geom,nion) 53 if (.not. value) call errquit('opening qmmgeometry',0, 54 & GEOM_ERR) 55 nwater = nion/3 56 else 57 nwater = 0 58 end if 59 60 61 if (nwater.gt.0) then 62 63* ***** Boundary condtion check **** 64 if (control_version().ne.4) then 65 value = geom_destroy(geom) 66 call errquit( 67 > 'Aperiodic boundary conditions must be used with Waterpsp',0, 68 & GEOM_ERR) 69 end if 70 71* ***** allocate waterpsp data structure ***** 72 value = BA_alloc_get(mt_dbl,(3*nion), 73 . 'rwater2',rwater2(2),rwater2(1)) 74 value = value.and. 75 > BA_alloc_get(mt_dbl,(3*nion), 76 > 'rwater1',rwater1(2),rwater1(1)) 77 value = value.and. 78 > BA_alloc_get(mt_dbl,(3*nion), 79 > 'rwater0',rwater0(2),rwater0(1)) 80 if (.not. value) call errquit('out of heap memory',0, MA_ERR) 81 82 do i=1,nion 83 value = value.and. 84 > geom_cent_get(geom,i,t,dbl_mb(rwater1(1)+(i-1)*3),q) 85 end do 86 call dcopy(3*nion,dbl_mb(rwater1(1)),1,dbl_mb(rwater2(1)),1) 87 88 call dcopy(3*nion,0.0d0,0,dbl_mb(rwater0(1)),1) 89 value = value.and.geom_vel_get(geom,dbl_mb(rwater0(1))) 90 value = value.and.geom_destroy(geom) 91 if (.not. value) call errquit('error reading qmmmgeometry', 0, 92 & GEOM_ERR) 93 call dcopy(3*nion,0.0d0,0,dbl_mb(rwater0(1)),1) 94 95* *** Init BLCJ SR Potential and LJparam **** 96 call LJparam_init(rtdb) 97 call BLCJ_SR_init(rtdb) 98 call BLCJ_LR_init(rtdb) 99 end if 100 101* **** set ekw_total,ekw_count **** 102 ekw_total = 0.0d0 103 ekw_count = 0 104 105 return 106 end 107 108 109* ********************************** 110* * * 111* * Waterpsp_write * 112* * * 113* ********************************** 114 115 subroutine Waterpsp_write(rtdb) 116 implicit none 117 integer rtdb 118 119***** water common block **** 120#include "waterpsp.fh" 121#include "bafdecls.fh" 122#include "geom.fh" 123 124* **** local variables **** 125 logical value 126 integer i,geom,nion 127 double precision rxyz(3),q 128 character*16 t 129 130 if (nwater.gt.0) then 131 value = geom_create(geom,'qmmmgeometry') 132 value = geom_rtdb_load(rtdb,geom,'qmmmgeometry') 133 value = geom_ncent(geom,nion) 134 do i=1,nion 135 value = geom_cent_get(geom,i,t,rxyz,q) 136 value = geom_cent_set(geom,i,t,dbl_mb(rwater1(1)+(i-1)*3),q) 137 end do 138 value = geom_vel_set(geom,dbl_mb(rwater0(1))) 139 140 value = geom_rtdb_delete(rtdb,'qmmmgeometry') 141 value = geom_rtdb_store(rtdb,geom,'qmmmgeometry') 142 value = geom_destroy(geom) 143 end if 144 145 return 146 end 147 148* ********************************** 149* * * 150* * Waterpsp_found * 151* * * 152* ********************************** 153 logical function Waterpsp_found() 154 implicit none 155 156***** water common block **** 157#include "waterpsp.fh" 158 159 logical value 160 161 value = .false. 162 if (nwater.gt.0) value = .true. 163 164 Waterpsp_found = value 165 return 166 end 167 168* ********************************** 169* * * 170* * Waterpsp_end * 171* * * 172* ********************************** 173 174 subroutine Waterpsp_end() 175 implicit none 176 177***** water common block **** 178#include "waterpsp.fh" 179#include "bafdecls.fh" 180 181 logical value 182 183 if (nwater.gt.0) then 184 value = BA_free_heap(rwater2(2)) 185 value = BA_free_heap(rwater1(2)) 186 value = BA_free_heap(rwater0(2)) 187 call LJparam_end() 188 call BLCJ_SR_end() 189 end if 190 191 return 192 end 193 194* ********************************** 195* * * 196* * Waterpsp_nwater * 197* * * 198* ********************************** 199 integer function Waterpsp_nwater() 200 implicit none 201 202* ***** water common block **** 203#include "waterpsp.fh" 204 205 Waterpsp_nwater = nwater 206 return 207 end 208 209* ********************************** 210* * * 211* * Waterpsp_ke * 212* * * 213* ********************************** 214 real*8 function Waterpsp_ke() 215 implicit none 216 217* ***** water common block **** 218#include "waterpsp.fh" 219 220 Waterpsp_ke = ekw 221 return 222 end 223 224* ********************************** 225* * * 226* * Waterpsp_Temperature * 227* * * 228* ********************************** 229 real*8 function Waterpsp_Temperature() 230 implicit none 231 232* ***** water common block **** 233#include "waterpsp.fh" 234 235 real*8 kb 236 parameter (kb=3.16679d-6) 237 238 Waterpsp_Temperature = 2.0d0*ekw_total/dble(ekw_count) 239 > /(9.0d0*nwater-6.0d0) 240 > /kb 241 return 242 end 243 244 245* ********************************** 246* * * 247* * Waterpsp_Print * 248* * * 249* ********************************** 250 251 subroutine Waterpsp_Print(unit) 252 implicit none 253 integer unit 254 255* ***** water common block **** 256#include "waterpsp.fh" 257#include "bafdecls.fh" 258 259* ***** local variables **** 260 integer taskid,MASTER 261 parameter (MASTER=0) 262 integer ii,k,i1,i2,i3 263 264 if (nwater.gt.0) then 265 call Parallel_taskid(taskid) 266 267 if (taskid.eq.MASTER) then 268 write(unit,1180) 269 do ii=1,nwater 270 i1 = 9*(ii-1) 271 i2 = i1+3 272 i3 = i2+3 273 write(unit,1190) ii,'O^',(dbl_mb(rwater1(1)+k-1+i1),k=1,3) 274 write(unit,1190) ii,'H^',(dbl_mb(rwater1(1)+k-1+i2),k=1,3) 275 write(unit,1190) ii,'H^',(dbl_mb(rwater1(1)+k-1+i3),k=1,3) 276 end do 277 end if 278 279 end if 280 return 281 282 1180 FORMAT(/' position of pseudopotential water:') 283 1190 FORMAT(5X, I4, A3 ,' (',3F11.5,' )') 284 end 285 286* ********************************** 287* * * 288* * Waterpsp_Print2 * 289* * * 290* ********************************** 291 292 subroutine Waterpsp_Print2(unit) 293 implicit none 294 integer unit 295 296* ***** water common block **** 297#include "waterpsp.fh" 298#include "bafdecls.fh" 299 300* ***** local variables **** 301 integer taskid,MASTER 302 parameter (MASTER=0) 303 integer ii,k,i1,i2,i3 304 305 if (nwater.gt.0) then 306 call Parallel_taskid(taskid) 307 308 if (taskid.eq.MASTER) then 309 write(unit,1180) 310 do ii=1,nwater 311 i1 = 9*(ii-1) 312 i2 = i1+3 313 i3 = i2+3 314 write(unit,1190) ii,'O^',(dbl_mb(rwater0(1)+k-1+i1),k=1,3) 315 write(unit,1190) ii,'H^',(dbl_mb(rwater0(1)+k-1+i2),k=1,3) 316 write(unit,1190) ii,'H^',(dbl_mb(rwater0(1)+k-1+i3),k=1,3) 317 end do 318 end if 319 320 end if 321 return 322 323 1180 FORMAT(/' velocity of pseudopotential water:') 324 1190 FORMAT(5X, I4, A3 ,' (',3F11.5,' )') 325 end 326 327 328* ********************************** 329* * * 330* * Waterpsp_PrintXYZ * 331* * * 332* ********************************** 333 334 subroutine Waterpsp_PrintXYZ(unit) 335 implicit none 336 integer unit 337 338* ***** water common block **** 339#include "waterpsp.fh" 340#include "bafdecls.fh" 341 342* ***** local variables **** 343 integer taskid,MASTER 344 parameter (MASTER=0) 345 integer ii,k,i1,i2,i3 346 347 if (nwater.gt.0) then 348 call Parallel_taskid(taskid) 349 350 if (taskid.eq.MASTER) then 351 do ii=1,nwater 352 i1 = 9*(ii-1) 353 i2 = i1+3 354 i3 = i2+3 355 write(unit,*) 'O ',(dbl_mb(rwater1(1)+k-1+i1)*0.529177d0,k=1,3) 356 write(unit,*) 'H ',(dbl_mb(rwater1(1)+k-1+i2)*0.529177d0,k=1,3) 357 write(unit,*) 'H ',(dbl_mb(rwater1(1)+k-1+i3)*0.529177d0,k=1,3) 358 end do 359 end if 360 361 end if 362 363 return 364 end 365 366* ********************************** 367* * * 368* * Waterpsp_rwater * 369* * * 370* ********************************** 371 real*8 function Waterpsp_rwater(i,j,ii) 372 implicit none 373 integer i,j,ii 374 375* ***** water common block **** 376#include "waterpsp.fh" 377#include "bafdecls.fh" 378 379 integer index 380 index = 9*(ii-1) + 3*(j-1) + (i-1) 381 Waterpsp_rwater = dbl_mb(rwater1(1)+index) 382 return 383 end 384 385* ********************************** 386* * * 387* * Waterpsp_vwater * 388* * * 389* ********************************** 390 real*8 function Waterpsp_vwater(i,j,ii) 391 implicit none 392 integer i,j,ii 393 394* ***** water common block **** 395#include "waterpsp.fh" 396#include "bafdecls.fh" 397 398 integer index 399 index = 9*(ii-1) + 3*(j-1) + (i-1) 400 Waterpsp_vwater = dbl_mb(rwater0(1)+index) 401 return 402 end 403 404 405* ********************************** 406* * * 407* * Waterpsp_shift * 408* * * 409* ********************************** 410 subroutine Waterpsp_shift() 411 implicit none 412 413* ***** water common block **** 414#include "waterpsp.fh" 415#include "bafdecls.fh" 416 417 call dcopy((9*nwater),dbl_mb(rwater1(1)),1,dbl_mb(rwater0(1)),1) 418 call dcopy((9*nwater),dbl_mb(rwater2(1)),1,dbl_mb(rwater1(1)),1) 419 420 return 421 end 422 423* ********************************** 424* * * 425* * Waterpsp_Generate_V * 426* * * 427* ********************************** 428 429 subroutine Waterpsp_Generate_V(n2ft3d,rgrid,Vwpsp) 430 implicit none 431 integer n2ft3d 432 real*8 rgrid(3,*) 433 real*8 Vwpsp(*) 434 435* ***** water common block **** 436#include "waterpsp.fh" 437#include "bafdecls.fh" 438 439* **** local variables **** 440 integer ii,i1,i2,i3 441 442 443 if (nwater.gt.0) then 444 do ii=1,nwater 445 i1 = 9*(ii-1) 446 i2 = i1+3 447 i3 = i2+3 448 call BLCJ_LR(dbl_mb(rwater1(1)+i1), 449 > dbl_mb(rwater1(1)+i2), 450 > dbl_mb(rwater1(1)+i3), 451 > n2ft3d,rgrid, 452 > Vwpsp) 453 call BLCJ_SR(dbl_mb(rwater1(1)+i1), 454 > dbl_mb(rwater1(1)+i2), 455 > dbl_mb(rwater1(1)+i3), 456 > n2ft3d,rgrid, 457 > Vwpsp) 458 end do 459 end if 460 461 return 462 end 463 464 465 466* ********************************** 467* * * 468* * Waterpsp_Energy_Water * 469* * * 470* ********************************** 471 472 subroutine Waterpsp_Energy_Water(ewater) 473 implicit none 474 real*8 ewater 475 476* ***** water common block **** 477#include "waterpsp.fh" 478#include "bafdecls.fh" 479 480* **** local variables **** 481 integer ii,jj 482 integer i1,i2,i3 483 integer j1,j2,j3 484 485* **** external functions **** 486 real*8 BLCJ_Intramolecular,BLCJ_Intermolecular 487 external BLCJ_Intramolecular,BLCJ_Intermolecular 488 489 490 ewater = 0.0d0 491 if (nwater.gt.0) then 492* **** calculate Intramolecular contributions **** 493 do ii=1,nwater 494 i1 = 9*(ii-1) 495 i2 = i1+3 496 i3 = i2+3 497 ewater = ewater 498 > + BLCJ_Intramolecular(dbl_mb(rwater1(1)+i1), 499 > dbl_mb(rwater1(1)+i2), 500 > dbl_mb(rwater1(1)+i3)) 501 502 end do 503 504* **** calculate Intermolecular contributions **** 505 do ii=1,nwater-1 506 i1 = 9*(ii-1) 507 i2 = i1+3 508 i3 = i2+3 509 do jj=ii+1,nwater 510 j1 = 9*(jj-1) 511 j2 = j1+3 512 j3 = j2+3 513 ewater = ewater 514 > + BLCJ_Intermolecular(dbl_mb(rwater1(1)+i1), 515 > dbl_mb(rwater1(1)+i2), 516 > dbl_mb(rwater1(1)+i3), 517 > dbl_mb(rwater1(1)+j1), 518 > dbl_mb(rwater1(1)+j2), 519 > dbl_mb(rwater1(1)+j3)) 520 end do 521 end do 522 523 end if 524 525 return 526 end 527 528* ********************************** 529* * * 530* * Waterpsp_Energy_ion * 531* * * 532* ********************************** 533 534 subroutine Waterpsp_Energy_ion(ewater) 535 implicit none 536 537 real*8 ewater 538 539* ***** water common block **** 540#include "waterpsp.fh" 541#include "bafdecls.fh" 542 543* **** local variables **** 544 integer ii,j,i1,i2,i3 545 real*8 zv,ec,er,ep 546 real*8 ei,si,ew,sw 547 real*8 rion(3) 548 549* **** external functions **** 550 integer ion_nion,ion_katm 551 real*8 psp_zv,ion_rion 552 external ion_nion,ion_katm 553 external psp_zv,ion_rion 554 555 ewater = 0.0d0 556 557 if (nwater.gt.0) then 558 call LJparam_water(ew,sw) 559 do j=1,ion_nion() 560 zv = psp_zv(ion_katm(j)) 561 rion(1) = ion_rion(1,j) 562 rion(2) = ion_rion(2,j) 563 rion(3) = ion_rion(3,j) 564 call LJparam_ion(j,ei,si) 565 do ii=1,nwater 566 i1 = 9*(ii-1) 567 i2 = i1+3 568 i3 = i2+3 569 ec = 0.0d0 570 call BLCJ_ion_Coulomb(dbl_mb(rwater1(1)+i1), 571 > dbl_mb(rwater1(1)+i2), 572 > dbl_mb(rwater1(1)+i3), 573 > rion,zv, 574 > ec) 575 ep = 0.0d0 576 call BLCJ_ion_Polarization(dbl_mb(rwater1(1)+i1), 577 > dbl_mb(rwater1(1)+i2), 578 > dbl_mb(rwater1(1)+i3), 579 > rion,zv, 580 > ep) 581 er = 0.0d0 582 call LJ_ion_Repulsion(dbl_mb(rwater1(1)+i1),ew,sw, 583 > rion,ei,si,er) 584 585 ewater = ewater + ec + ep + er 586 end do 587 end do 588 589 end if 590 591 return 592 end 593 594* ********************************** 595* * * 596* * Waterpsp_Fion * 597* * * 598* ********************************** 599 subroutine Waterpsp_Fion(fion) 600 implicit none 601 real*8 fion(3,*) 602 603* ***** water common block **** 604#include "waterpsp.fh" 605#include "bafdecls.fh" 606 607* **** local variables **** 608 integer ii,j,i1,i2,i3 609 real*8 fc(3),fp(3),fl(3) 610 real*8 zv,ei,si,ew,sw 611 real*8 rion(3) 612 613 614* **** external functions **** 615 integer ion_katm,ion_nion 616 real*8 psp_zv,ion_rion 617 external ion_katm,ion_nion 618 external psp_zv,ion_rion 619 620 621* **** calculate the force on QM ions **** 622 call LJparam_water(ew,sw) 623 do j=1,ion_nion() 624 zv = psp_zv(ion_katm(j)) 625 rion(1) = ion_rion(1,j) 626 rion(2) = ion_rion(2,j) 627 rion(3) = ion_rion(3,j) 628 call LJparam_ion(j,ei,si) 629 630 do ii=1,nwater 631 i1 = 9*(ii-1) 632 i2 = i1+3 633 i3 = i2+3 634 635 call dcopy(3,0.0d0,0,fc,1) 636 call dcopy(3,0.0d0,0,fp,1) 637 call dcopy(3,0.0d0,0,fl,1) 638 call BLCJ_ion_Coulomb_Fion(dbl_mb(rwater1(1) + i1), 639 > dbl_mb(rwater1(1) + i2), 640 > dbl_mb(rwater1(1) + i3), 641 > rion,zv,fc) 642 call BLCJ_ion_Polarization_Fion(dbl_mb(rwater1(1)+i1), 643 > dbl_mb(rwater1(1)+i2), 644 > dbl_mb(rwater1(1)+i3), 645 > rion,zv,fp) 646 call LJ_ion_Repulsion_Force(dbl_mb(rwater1(1)+i1),ew,sw, 647 > rion,ei,si,fl) 648 649 fion(1,j) = fion(1,j) + fc(1) + fp(1) + fl(1) 650 fion(2,j) = fion(2,j) + fc(2) + fp(2) + fl(2) 651 fion(3,j) = fion(3,j) + fc(3) + fp(3) + fl(3) 652 end do 653 654 end do 655 return 656 end 657 658 659* ********************************** 660* * * 661* * Waterpsp_Fwater * 662* * * 663* ********************************** 664 subroutine Waterpsp_Fwater(n2ft3d,rgrid,rho,dv,fwater) 665 implicit none 666 integer n2ft3d 667 real*8 rgrid(3,*) 668 real*8 rho(*) 669 real*8 dv 670 real*8 fwater(3,3,*) 671 672 673* ***** water common block **** 674#include "waterpsp.fh" 675#include "bafdecls.fh" 676 677* **** local variables **** 678 integer j,ii,i1,i2,i3 679 integer jj,j1,j2,j3 680 real*8 ew,sw,ei,si,zv,rion(3) 681 real*8 foww(3),f1ww(3),f2ww(3) 682 real*8 foiw(3),f1iw(3),f2iw(3) 683 real*8 foew(3),f1ew(3),f2ew(3) 684 real*8 foc(3),f1c(3),f2c(3),fop(3) 685 real*8 fol(3),f1l(3),f2l(3) 686 687* **** external functions **** 688 integer ion_katm,ion_nion 689 real*8 psp_zv,ion_rion 690 external ion_katm,ion_nion 691 external psp_zv,ion_rion 692 693 call LJparam_water(ew,sw) 694 695* **** calculation the force on waters **** 696 do ii=1,nwater 697 i1 = 9*(ii-1) 698 i2 = i1+3 699 i3 = i2+3 700 701* **** Intra water-water interaction **** 702 call dcopy(3,0.0d0,0,foww,1) 703 call dcopy(3,0.0d0,0,f1ww,1) 704 call dcopy(3,0.0d0,0,f2ww,1) 705 call BLCJ_Intramolecular_Fwater(dbl_mb(rwater1(1)+i1), 706 > dbl_mb(rwater1(1)+i2), 707 > dbl_mb(rwater1(1)+i3), 708 > foww,f1ww,f2ww) 709 710* **** ion-water interaction **** 711 call dcopy(3,0.0d0,0,foiw,1) 712 call dcopy(3,0.0d0,0,f1iw,1) 713 call dcopy(3,0.0d0,0,f2iw,1) 714 do j=1,ion_nion() 715 zv = psp_zv(ion_katm(j)) 716 rion(1) = ion_rion(1,j) 717 rion(2) = ion_rion(2,j) 718 rion(3) = ion_rion(3,j) 719 call LJparam_ion(j,ei,si) 720 721 call dcopy(3,0.0d0,0,foc,1) 722 call dcopy(3,0.0d0,0,f1c,1) 723 call dcopy(3,0.0d0,0,f2c,1) 724 call dcopy(3,0.0d0,0,fop,1) 725 call dcopy(3,0.0d0,0,fol,1) 726 call BLCJ_ion_Coulomb_Fwater(dbl_mb(rwater1(1)+i1), 727 > dbl_mb(rwater1(1)+i2), 728 > dbl_mb(rwater1(1)+i3), 729 > rion,zv,foc,f1c,f2c) 730 call BLCJ_ion_Polarization_Fwater(dbl_mb(rwater1(1)+i1), 731 > dbl_mb(rwater1(1)+i2), 732 > dbl_mb(rwater1(1)+i3), 733 > rion,zv,fop) 734 call LJ_ion_Repulsion_Force(rion,ei,si, 735 > dbl_mb(rwater1(1)+i1),ew,sw, 736 > fol) 737 738 foiw(1) = foiw(1) + foc(1) + fop(1) + fol(1) 739 foiw(2) = foiw(2) + foc(2) + fop(2) + fol(2) 740 foiw(3) = foiw(3) + foc(3) + fop(3) + fol(3) 741 742 f1iw(1) = f1iw(1) + f1c(1) 743 f1iw(2) = f1iw(2) + f1c(2) 744 f1iw(3) = f1iw(3) + f1c(3) 745 746 f2iw(1) = f2iw(1) + f2c(1) 747 f2iw(2) = f2iw(2) + f2c(2) 748 f2iw(3) = f2iw(3) + f2c(3) 749 750 end do 751 752* **** elc-water interaction **** 753 call dcopy(3,0.0d0,0,foc,1) 754 call dcopy(3,0.0d0,0,f1c,1) 755 call dcopy(3,0.0d0,0,f2c,1) 756 call dcopy(3,0.0d0,0,fol,1) 757 call dcopy(3,0.0d0,0,f1l,1) 758 call dcopy(3,0.0d0,0,f2l,1) 759 call BLCJ_LR_Fwater(dbl_mb(rwater1(1)+i1), 760 > dbl_mb(rwater1(1)+i2), 761 > dbl_mb(rwater1(1)+i3), 762 > n2ft3d,rgrid,rho,dv, 763 > foc,f1c,f2c) 764 call BLCJ_SR_Fwater(dbl_mb(rwater1(1)+i1), 765 > dbl_mb(rwater1(1)+i2), 766 > dbl_mb(rwater1(1)+i3), 767 > n2ft3d,rgrid,rho,dv, 768 > fol,f1l,f2l) 769 770 foew(1) = foc(1) + fol(1) 771 foew(2) = foc(2) + fol(2) 772 foew(3) = foc(3) + fol(3) 773 774 f1ew(1) = f1c(1) + f1l(1) 775 f1ew(2) = f1c(2) + f1l(2) 776 f1ew(3) = f1c(3) + f1l(3) 777 778 f2ew(1) = f2c(1) + f2l(1) 779 f2ew(2) = f2c(2) + f2l(2) 780 f2ew(3) = f2c(3) + f2l(3) 781 782* **** sum up contributions **** 783 fwater(1,1,ii) = foww(1) + foiw(1) + foew(1) 784 fwater(2,1,ii) = foww(2) + foiw(2) + foew(2) 785 fwater(3,1,ii) = foww(3) + foiw(3) + foew(3) 786 787 fwater(1,2,ii) = f1ww(1) + f1iw(1) + f1ew(1) 788 fwater(2,2,ii) = f1ww(2) + f1iw(2) + f1ew(2) 789 fwater(3,2,ii) = f1ww(3) + f1iw(3) + f1ew(3) 790 791 fwater(1,3,ii) = f2ww(1) + f2iw(1) + f2ew(1) 792 fwater(2,3,ii) = f2ww(2) + f2iw(2) + f2ew(2) 793 fwater(3,3,ii) = f2ww(3) + f2iw(3) + f2ew(3) 794 795 end do 796 797* **** Inter water-water interaction **** 798 do ii=1,nwater 799 i1 = 9*(ii-1) 800 i2 = i1+3 801 i3 = i2+3 802 do jj=ii+1,nwater 803 j1 = 9*(jj-1) 804 j2 = j1+3 805 j3 = j2+3 806 call BLCJ_Intermolecular_Fwater(dbl_mb(rwater1(1)+i1), 807 > dbl_mb(rwater1(1)+i2), 808 > dbl_mb(rwater1(1)+i3), 809 > fwater(1,1,ii),fwater(1,2,ii),fwater(1,3,ii), 810 > dbl_mb(rwater1(1)+j1), 811 > dbl_mb(rwater1(1)+j2), 812 > dbl_mb(rwater1(1)+j3), 813 > fwater(1,1,jj),fwater(1,2,jj),fwater(1,3,jj)) 814 end do 815 end do 816 817 return 818 end 819 820 821 822* ********************************** 823* * * 824* * Waterpsp_Update * 825* * * 826* ********************************** 827 subroutine Waterpsp_Update(algorithm, 828 > n2ft3d,rgrid,rho,dv, 829 > dt,fion,deltawater) 830 implicit none 831#include "errquit.fh" 832 integer algorithm 833 integer n2ft3d 834 real*8 rgrid(3,*) 835 real*8 rho(*) 836 real*8 dv 837 real*8 dt 838 real*8 fion(3,*) 839 real*8 deltawater 840 841* ***** water common block **** 842#include "waterpsp.fh" 843#include "bafdecls.fh" 844 845* **** local variables **** 846 logical value 847 integer fwater(2) 848 849 850* **** allocate fwater force *** 851 value = BA_push_get(mt_dbl,(9*nwater), 852 > 'fwater',fwater(2),fwater(1)) 853 if (.not.value) 854 > call errquit('Waterpsp_SD_Update:error push stack',0, MA_ERR) 855 856 857* **** Get Forces **** 858 call Waterpsp_Fion(fion) 859 call Waterpsp_Fwater(n2ft3d,rgrid,rho,dv,dbl_mb(fwater(1))) 860 deltawater = 0.0d0 861 862 863* **** steepest descent update of water **** 864 if (algorithm.eq.0) then 865 call Waterpsp_SD_subUpdate(nwater,dbl_mb(rwater2(1)), 866 > dbl_mb(rwater1(1)), 867 > dbl_mb(fwater(1)), 868 > dt, 869 > deltawater) 870 871* **** Newton step update of water **** 872 else if (algorithm.eq.1) then 873 call Waterpsp_Newton_subUpdate(nwater,dbl_mb(rwater2(1)), 874 > dbl_mb(rwater1(1)), 875 > dbl_mb(rwater0(1)), 876 > dbl_mb(fwater(1)), 877 > dt,ekw) 878 ekw_total = ekw_total + ekw 879 ekw_count = ekw_count + 1 880 881* **** Verlet step update of water **** 882 else if (algorithm.eq.2) then 883 call Waterpsp_Verlet_subUpdate(nwater,dbl_mb(rwater2(1)), 884 > dbl_mb(rwater1(1)), 885 > dbl_mb(rwater0(1)), 886 > dbl_mb(fwater(1)), 887 > dt,ekw) 888 ekw_total = ekw_total + ekw 889 ekw_count = ekw_count + 1 890 end if 891 892 893* **** deallocate fwater force *** 894 value = BA_pop_stack(fwater(2)) 895 if (.not.value) 896 > call errquit('Waterpsp_SD_Update:error pop stack',0, MA_ERR) 897 return 898 end 899 900 901 902* ********************************** 903* * * 904* * Waterpsp_SD_subUpdate * 905* * * 906* ********************************** 907 subroutine Waterpsp_SD_subUpdate(nwater,rwater2,rwater1,fwater, 908 > dt,deltawater) 909 implicit none 910 integer nwater 911 real*8 rwater2(3,3,*) 912 real*8 rwater1(3,3,*) 913 real*8 fwater(3,3,*) 914 real*8 dt 915 real*8 deltawater 916 917* **** local variables *** 918 integer j,ii 919 real*8 dtO,dtH,sum 920 921* **** steepest descent update of water **** 922 dtO = dt/dsqrt(16.0d0*1822.89d0) 923 dtH = dt/dsqrt( 1.0d0*1822.89d0) 924 do ii=1,nwater 925 rwater2(1,1,ii) = rwater1(1,1,ii) + dtO*fwater(1,1,ii) 926 rwater2(2,1,ii) = rwater1(2,1,ii) + dtO*fwater(2,1,ii) 927 rwater2(3,1,ii) = rwater1(3,1,ii) + dtO*fwater(3,1,ii) 928 929 rwater2(1,2,ii) = rwater1(1,2,ii) + dtH*fwater(1,2,ii) 930 rwater2(2,2,ii) = rwater1(2,2,ii) + dtH*fwater(2,2,ii) 931 rwater2(3,2,ii) = rwater1(3,2,ii) + dtH*fwater(3,2,ii) 932 933 rwater2(1,3,ii) = rwater1(1,3,ii) + dtH*fwater(1,3,ii) 934 rwater2(2,3,ii) = rwater1(2,3,ii) + dtH*fwater(2,3,ii) 935 rwater2(3,3,ii) = rwater1(3,3,ii) + dtH*fwater(3,3,ii) 936 end do 937 938* **** find maximum force **** 939 deltawater = 0.0d0 940 do ii=1,nwater 941 do j=1,3 942 sum = dsqrt(fwater(1,j,ii)**2 943 > +fwater(2,j,ii)**2 944 > +fwater(3,j,ii)**2) 945 if (sum.gt.deltawater) deltawater = sum 946 end do 947 end do 948 949 return 950 end 951 952* ********************************** 953* * * 954* * Waterpsp_Newton_subUpdate * 955* * * 956* ********************************** 957 subroutine Waterpsp_Newton_subUpdate(nwater,rwater2,rwater1, 958 > vwater,fwater,dt,ekw) 959 implicit none 960 integer nwater 961 real*8 rwater2(3,3,*) 962 real*8 rwater1(3,3,*) 963 real*8 vwater(3,3,*) 964 real*8 fwater(3,3,*) 965 real*8 dt,ekw 966 967* **** local variables **** 968 integer ii 969 real*8 dtO,dtH,MO,MH 970 971* **** steepest descent update of water **** 972 dtO = 0.5d0*dt*dt/(16.0d0*1822.89d0) 973 dtH = 0.5d0*dt*dt/( 1.0d0*1822.89d0) 974 do ii=1,nwater 975 rwater2(1,1,ii) = rwater1(1,1,ii) + dt *vwater(1,1,ii) 976 > + dtO*fwater(1,1,ii) 977 rwater2(2,1,ii) = rwater1(2,1,ii) + dt *vwater(2,1,ii) 978 > + dtO*fwater(2,1,ii) 979 rwater2(3,1,ii) = rwater1(3,1,ii) + dt *vwater(3,1,ii) 980 > + dtO*fwater(3,1,ii) 981 982 rwater2(1,2,ii) = rwater1(1,2,ii) + dt *vwater(1,2,ii) 983 > + dtH*fwater(1,2,ii) 984 rwater2(2,2,ii) = rwater1(2,2,ii) + dt *vwater(2,2,ii) 985 > + dtH*fwater(2,2,ii) 986 rwater2(3,2,ii) = rwater1(3,2,ii) + dt *vwater(3,2,ii) 987 > + dtH*fwater(3,2,ii) 988 989 rwater2(1,3,ii) = rwater1(1,3,ii) + dt *vwater(1,3,ii) 990 > + dtH*fwater(1,3,ii) 991 rwater2(2,3,ii) = rwater1(2,3,ii) + dt *vwater(2,3,ii) 992 > + dtH*fwater(2,3,ii) 993 rwater2(3,3,ii) = rwater1(3,3,ii) + dt *vwater(3,3,ii) 994 > + dtH*fwater(3,3,ii) 995 end do 996 997* ***** find kinetic energy **** 998 ekw = 0.0d0 999 MO = (16.0d0*1822.89d0) 1000 MH = ( 1.0d0*1822.89d0) 1001 do ii=1,nwater 1002 ekw = ekw + MO*(vwater(1,1,ii)**2 1003 > +vwater(2,1,ii)**2 1004 > +vwater(3,1,ii)**2) 1005 > + MH*(vwater(1,2,ii)**2 1006 > +vwater(2,2,ii)**2 1007 > +vwater(3,2,ii)**2) 1008 > + MH*(vwater(1,3,ii)**2 1009 > +vwater(2,3,ii)**2 1010 > +vwater(3,3,ii)**2) 1011 end do 1012 ekw = 0.5d0*ekw 1013 1014 return 1015 end 1016 1017 1018* ********************************** 1019* * * 1020* * Waterpsp_Verlet_subUpdate * 1021* * * 1022* ********************************** 1023 subroutine Waterpsp_Verlet_subUpdate(nwater, 1024 > rwater2,rwater1,rwater0,fwater, 1025 > dt,ekw) 1026 implicit none 1027 integer nwater 1028 real*8 rwater2(3,3,*) 1029 real*8 rwater1(3,3,*) 1030 real*8 rwater0(3,3,*) 1031 real*8 fwater(3,3,*) 1032 real*8 dt,ekw 1033 1034* **** local variables **** 1035 integer i,j,ii 1036 real*8 h,dtO,dtH,MO,MH 1037 1038* **** steepest descent update of water **** 1039 dtO = dt*dt/(16.0d0*1822.89d0) 1040 dtH = dt*dt/( 1.0d0*1822.89d0) 1041 do ii=1,nwater 1042 rwater2(1,1,ii) = 2*rwater1(1,1,ii) 1043 > - rwater0(1,1,ii) + dtO*fwater(1,1,ii) 1044 rwater2(2,1,ii) = 2*rwater1(2,1,ii) 1045 > - rwater0(2,1,ii) + dtO*fwater(2,1,ii) 1046 rwater2(3,1,ii) = 2*rwater1(3,1,ii) 1047 > - rwater0(3,1,ii) + dtO*fwater(3,1,ii) 1048 1049 rwater2(1,2,ii) = 2*rwater1(1,2,ii) 1050 > - rwater0(1,2,ii) + dtH*fwater(1,2,ii) 1051 rwater2(2,2,ii) = 2*rwater1(2,2,ii) 1052 > - rwater0(2,2,ii) + dtH*fwater(2,2,ii) 1053 rwater2(3,2,ii) = 2*rwater1(3,2,ii) 1054 > - rwater0(3,2,ii) + dtH*fwater(3,2,ii) 1055 1056 rwater2(1,3,ii) = 2*rwater1(1,3,ii) 1057 > - rwater0(1,3,ii) + dtH*fwater(1,3,ii) 1058 rwater2(2,3,ii) = 2*rwater1(2,3,ii) 1059 > - rwater0(2,3,ii) + dtH*fwater(2,3,ii) 1060 rwater2(3,3,ii) = 2*rwater1(3,3,ii) 1061 > - rwater0(3,3,ii) + dtH*fwater(3,3,ii) 1062 end do 1063 1064* **** make rwater0 the velocity - note that velocity is deleted after **** 1065* **** Waterpsp_shift call **** 1066 h = 1.0d0/(2.0d0*dt) 1067 do ii=1,nwater 1068 do j=1,3 1069 do i=1,3 1070 rwater0(i,j,ii) = h*(rwater2(i,j,ii) - rwater0(i,j,ii)) 1071 end do 1072 end do 1073 end do 1074 1075* ***** find kinetic energy **** 1076 ekw = 0.0d0 1077 MO = (16.0d0*1822.89d0) 1078 MH = ( 1.0d0*1822.89d0) 1079 do ii=1,nwater 1080 ekw = ekw + MO*(rwater0(1,1,ii)**2 1081 > +rwater0(2,1,ii)**2 1082 > +rwater0(3,1,ii)**2) 1083 > + MH*(rwater0(1,2,ii)**2 1084 > +rwater0(2,2,ii)**2 1085 > +rwater0(3,2,ii)**2) 1086 > + MH*(rwater0(1,3,ii)**2 1087 > +rwater0(2,3,ii)**2 1088 > +rwater0(3,3,ii)**2) 1089 end do 1090 ekw = 0.5d0*ekw 1091 1092 return 1093 end 1094c $Id$ 1095