1 subroutine smd_vlist_init_system() 2 implicit none 3#include "errquit.fh" 4#include "inp.fh" 5#include "mafdecls.fh" 6#include "rtdb.fh" 7#include "util.fh" 8#include "global.fh" 9c 10 character*32 sp_vlist,sp_exlist,sp_coords 11 character*32 tag,pname 12 logical result 13 logical oexlist 14 15 pname = "smd_vlist_init_system" 16c 17 tag = "coordinates" 18 call smd_system_get_component(sp_coords,tag,result) 19 if(.not.result) 20 > call errquit( 21 > pname//'no component '//tag,0,0) 22 23 oexlist = .true. 24 tag = "excl_list" 25 call smd_system_get_component(sp_exlist,tag,result) 26 if(.not.result) oexlist = .false. 27 28 tag = "verlet_list" 29 call smd_system_get_component(sp_vlist,tag,result) 30 if(.not.result) 31 > call errquit( 32 > pname//'no component '//tag,0,0) 33 34 35 call smd_vlist_init(oexlist,sp_vlist,sp_exlist,sp_coords) 36c 37 return 38 end 39 40 subroutine smd_vlist_init(oexlist,sp_vlist,sp_exlist,sp_coord) 41 implicit none 42#include "errquit.fh" 43#include "inp.fh" 44#include "mafdecls.fh" 45#include "rtdb.fh" 46#include "util.fh" 47#include "global.fh" 48c 49 character*(*) sp_vlist 50 character*(*) sp_exlist 51 character*(*) sp_coord 52 logical oexlist 53c 54 character*32 pname 55 character*80 tag 56 integer np,nl 57 integer i_p,i_l,i_c 58 integer h_l 59 integer i_list,i_clist 60 integer i 61 integer i_xp,i_xl 62 integer i_cl,h_cl 63 integer i_ct,h_ct 64 double precision rc2 65 integer nlb 66 logical result 67c 68 pname = "smd_vlist_init" 69c 70c write(*,*) "in "//pname 71c 72c get coordinate information 73c -------------------------- 74 tag = "coords" 75 call smd_data_get_index(sp_coord,tag,i_c,result) 76 if(.not. result) 77 > call errquit( 78 > pname//'error getting index for'//tag,0, RTDB_ERR) 79 80 call smd_data_get_size(sp_coord,tag,np,result) 81 if(.not. result) 82 > call errquit( 83 > pname//'error getting size for'//tag,0, RTDB_ERR) 84 np = np/3 85c 86c get excluded list information 87c ----------------------------- 88 if(oexlist) then 89 tag = "exlist:pointer" 90 call smd_data_get_index(sp_exlist,tag,i_xp,result) 91 if(.not. result) 92 > call errquit( 93 > pname//'error getting index for'//tag,0, RTDB_ERR) 94 tag = "exlist:list" 95 call smd_data_get_index(sp_exlist,tag,i_xl,result) 96 if(.not. result) 97 > call errquit( 98 > pname//'error getting index for'//tag,0, RTDB_ERR) 99 100 end if 101c 102c gestimate the size of pair list 103c ------------------------------ 104 nl = min( 7*np*200, ma_inquire_avail(MT_INT)) 105 nl = nl/7 106c 107c create pointer memory 108c --------------------- 109 call smd_namespace_create(sp_vlist) 110 tag = "vlist:pointer" 111 call smd_data_create(sp_vlist,"vlist:pointer",np,MT_INT) 112 call smd_data_get_index(sp_vlist,tag,i_p,result) 113 if(.not. result) 114 > call errquit( 115 > pname//'error getting index for'//tag,0, RTDB_ERR) 116 117c 118c create temporary scratch array for list since 119c we do not know the size yet 120c --------------------------------------------- 121 if(.not.ma_push_get(mt_int,nl,'tmp l',h_l,i_l)) 122 + call errquit(pname//'Failed to allocate memory for tmp l', 123 + nl, MA_ERR) 124 125 if(.not.ma_push_get(mt_dbl,3*nl,'tmp cl',h_cl,i_cl)) 126 + call errquit(pname//'Failed to allocate memory for tmp l', 127 + nl, MA_ERR) 128 129 if(.not.ma_push_get(mt_dbl,3*np,'tmp',h_ct,i_ct)) 130 + call errquit(pname//'Failed to allocate memory for tmp l', 131 + nl, MA_ERR) 132 133 134 call smd_cutoff_get_rcut_verlet(rc2) 135 rc2=rc2*rc2 136 137 if(oexlist) then 138 call smd_vlist_set(np, 139 + nl, 140 + rc2, 141 + dbl_mb(i_c), 142 + int_mb(i_xp), 143 + int_mb(i_xl), 144 + int_mb(i_p), 145 + int_mb(i_l), 146 + dbl_mb(i_cl), 147 + dbl_mb(i_ct), 148 + result) 149 150 else 151 152 call smd_vlist_set1(np, 153 + nl, 154 + rc2, 155 + dbl_mb(i_c), 156 + int_mb(i_p), 157 + int_mb(i_l), 158 + dbl_mb(i_cl), 159 + dbl_mb(i_ct), 160 + result) 161 162 163 end if 164c 165c create list memory 166c nl now contains the actual size 167c we will buffer it though to allow for possible expansion 168c -------------------------------------------------------- 169 nlb = 500 170 tag = "vlist:list" 171 call smd_data_create(sp_vlist,tag,nl+nlb,MT_INT) 172 call smd_data_get_index(sp_vlist,tag,i_list,result) 173 if(.not. result) 174 > call errquit( 175 > pname//'error getting index for'//tag,0, RTDB_ERR) 176 177 178 do i=1,nl 179 int_mb(i_list+i-1) = int_mb(i_l+i-1) 180 end do 181 182 tag = "vlist:distances" 183 call smd_data_create(sp_vlist,tag,3*(nl+nlb),MT_DBL) 184 call smd_data_get_index(sp_vlist,tag,i_clist,result) 185 if(.not. result) 186 > call errquit( 187 > pname//'error getting index for'//tag,0, RTDB_ERR) 188 189 tag = "vlist:displacement" 190 call smd_data_create(sp_vlist,tag,3*np,MT_DBL) 191 192 do i=1,3*nl 193 dbl_mb(i_clist+i-1) = dbl_mb(i_cl+i-1) 194 end do 195 196 if(.not.ma_pop_stack(h_ct)) 197 & call errquit(pname//'Failed to deallocate stack h_l',nl, 198 & MA_ERR) 199 200 201 if(.not.ma_pop_stack(h_cl)) 202 & call errquit(pname//'Failed to deallocate stack h_l',nl, 203 & MA_ERR) 204 205 206 if(.not.ma_pop_stack(h_l)) 207 & call errquit(pname//'Failed to deallocate stack h_l',nl, 208 & MA_ERR) 209 210 211 return 212 end 213c 214 subroutine smd_vlist_update(olist,ocoord) 215 implicit none 216#include "errquit.fh" 217#include "inp.fh" 218#include "mafdecls.fh" 219#include "rtdb.fh" 220#include "util.fh" 221#include "global.fh" 222c 223 logical olist 224 logical ocoord 225c 226 character*32 sp_vlist 227 character*32 sp_exlist 228 character*32 sp_coord 229c 230 character*32 pname 231 character*80 tag 232 integer np,nl 233 integer i_p,i_l,i_c 234 integer h_l 235 integer i_list,i_clist 236 integer i 237 integer i_xp,i_xl 238 double precision rc2 239 integer nlb 240 logical result 241 integer i_ct,h_ct 242 logical oexlist 243c 244 pname = "smd_vlist_update" 245c 246c write(*,*) "in "//pname 247c 248 if((.not.olist).and.(.not.ocoord)) return 249c 250c get components 251c -------------- 252 tag = "coordinates" 253 call smd_system_get_component(sp_coord,tag,result) 254 if(.not.result) 255 > call errquit( 256 > pname//'no component '//tag,0,0) 257 258 oexlist = .true. 259 tag = "excl_list" 260 call smd_system_get_component(sp_exlist,tag,result) 261 if(.not.result) oexlist =.false. 262 263 tag = "verlet_list" 264 call smd_system_get_component(sp_vlist,tag,result) 265 if(.not.result) 266 > call errquit( 267 > pname//'no component '//tag,0,0) 268 269c 270 if(.not.olist) then 271 call smd_vlist_update_coord(sp_vlist,sp_coord) 272 goto 200 273 end if 274c 275c get coordinate information 276c -------------------------- 277 tag = "coords" 278 call smd_data_get_index(sp_coord,tag,i_c,result) 279 if(.not. result) 280 > call errquit( 281 > pname//'error getting index for'//tag,0, RTDB_ERR) 282 283 call smd_data_get_size(sp_coord,tag,np,result) 284 if(.not. result) 285 > call errquit( 286 > pname//'error getting size for'//tag,0, RTDB_ERR) 287 np = np/3 288c 289c get excluded list information 290c ----------------------------- 291 if(oexlist) then 292 tag = "exlist:pointer" 293 call smd_data_get_index(sp_exlist,tag,i_xp,result) 294 if(.not. result) 295 > call errquit( 296 > pname//'error getting index for'//tag,0, RTDB_ERR) 297 tag = "exlist:list" 298 call smd_data_get_index(sp_exlist,tag,i_xl,result) 299 if(.not. result) 300 > call errquit( 301 > pname//'error getting index for'//tag,0, RTDB_ERR) 302 303 end if 304c 305c get verlet list information 306c --------------------------- 307 tag = "vlist:pointer" 308 call smd_data_get_index(sp_vlist,tag,i_p,result) 309 if(.not. result) 310 > call errquit( 311 > pname//'error getting index for'//tag,0, RTDB_ERR) 312 313 tag = "vlist:list" 314 call smd_data_get_index(sp_vlist,tag,i_list,result) 315 if(.not. result) 316 > call errquit( 317 > pname//'error getting index for'//tag,0, RTDB_ERR) 318 call smd_data_get_size(sp_vlist,tag,nl,result) 319 if(.not. result) 320 > call errquit( 321 > pname//'error getting size for'//tag,0, RTDB_ERR) 322 323 tag = "vlist:distances" 324 call smd_data_get_index(sp_vlist,tag,i_clist,result) 325 if(.not. result) 326 > call errquit( 327 > pname//'error getting index for'//tag,0, RTDB_ERR) 328c 329c create temporary scratch array for list 330c --------------------------------------------- 331 332 if(.not.ma_push_get(mt_dbl,3*np,'tmp',h_ct,i_ct)) 333 + call errquit(pname//'Failed to allocate memory for tmp l', 334 + nl, MA_ERR) 335 336 337 call smd_cutoff_get_rcut_verlet(rc2) 338 rc2=rc2*rc2 339 340 if(oexlist) then 341 call smd_vlist_set(np, 342 + nl, 343 + rc2, 344 + dbl_mb(i_c), 345 + int_mb(i_xp), 346 + int_mb(i_xl), 347 + int_mb(i_p), 348 + int_mb(i_list), 349 + dbl_mb(i_clist), 350 + dbl_mb(i_ct), 351 + result) 352 353 else 354 355 call smd_vlist_set1(np, 356 + nl, 357 + rc2, 358 + dbl_mb(i_c), 359 + int_mb(i_p), 360 + int_mb(i_list), 361 + dbl_mb(i_clist), 362 + dbl_mb(i_ct), 363 + result) 364 365 366 end if 367 368 if(.not.ma_pop_stack(h_ct)) 369 & call errquit(pname//'Failed to deallocate stack h_l',nl, 370 & MA_ERR) 371 372 373200 continue 374 return 375 end 376c 377 subroutine smd_vlist_update_coord(sp_vlist,sp_coord) 378 implicit none 379#include "errquit.fh" 380#include "inp.fh" 381#include "mafdecls.fh" 382#include "rtdb.fh" 383#include "util.fh" 384#include "global.fh" 385c 386 character*(*) sp_vlist 387 character*(*) sp_coord 388c 389 character*32 pname 390 character*80 tag 391 integer np,nl 392 integer i_p,i_l,i_c 393 integer i_list,i_clist 394 integer i_cl 395 logical result 396c 397 pname = "smd_vlist_init" 398c 399c write(*,*) "in "//pname 400c 401c get coordinate information 402c -------------------------- 403 tag = "coords" 404 call smd_data_get_index(sp_coord,tag,i_c,result) 405 if(.not. result) 406 > call errquit( 407 > pname//'error getting index for'//tag,0, RTDB_ERR) 408 409 call smd_data_get_size(sp_coord,tag,np,result) 410 if(.not. result) 411 > call errquit( 412 > pname//'error getting size for'//tag,0, RTDB_ERR) 413 np = np/3 414c 415c get list information 416c -------------------- 417 tag = "vlist:pointer" 418 call smd_data_get_index(sp_vlist,tag,i_p,result) 419 if(.not. result) 420 > call errquit( 421 > pname//'error getting index for'//tag,0, RTDB_ERR) 422 423 tag = "vlist:list" 424 call smd_data_get_index(sp_vlist,tag,i_list,result) 425 if(.not. result) 426 > call errquit( 427 > pname//'error getting index for'//tag,0, RTDB_ERR) 428 429 430 tag = "vlist:distances" 431 call smd_data_get_index(sp_vlist,tag,i_clist,result) 432 if(.not. result) 433 > call errquit( 434 > pname//'error getting index for'//tag,0, RTDB_ERR) 435 436 nl = int_mb(i_p+np-1) - int_mb(i_p) 437 call smd_vlist_update_coord0(np,nl, 438 + dbl_mb(i_c), 439 + int_mb(i_p), 440 + int_mb(i_l), 441 + dbl_mb(i_cl)) 442 443 call smd_lat_rebox(nl,dbl_mb(i_clist)) 444 445 return 446 end 447c 448 subroutine smd_vlist_update_coord0(np,nl,c,point, 449 > list,cl) 450c 451c update coordinates in verlet pair list 452c point(i) is an index to list() array 453c that contains all the pairs of atom i 454c In other words all the atoms paired with atom i 455c are contained in list(point(i)),...,list(point(i+1)-1) 456c np [in] number of atoms (which is also size of pointer array) 457c nl [in] list size 458c c [in] coordinates 459c point [in] verlet pointer 460c list [in] verlet list 461c cl [out] list of vectors (ri-rj) 462 implicit none 463#include "errquit.fh" 464 integer np,nl 465 double precision c(np,3) 466 integer point(np) 467 integer list(nl) 468 double precision cl(nl,3) 469c 470 integer i,j 471 integer nlist 472 integer jnab,jbeg,jend 473 474 character*30 pname 475 476 pname = "smd_vlist_update_coord0" 477 478 nlist=0 479 do i=1,np-1 480 jbeg=point(i) 481 jend=point(i+1)-1 482 483 if(jbeg.le.jend)then 484 485 do jnab=jbeg,jend 486 487 j=list(jnab) 488 489 nlist = nlist + 1 490 if(nlist.gt.nl) 491 > call errquit( 492 > pname//'out of bounds',0, RTDB_ERR) 493 494 cl(nlist,1)=c(i,1)-c(j,1) 495 cl(nlist,2)=c(i,2)-c(j,2) 496 cl(nlist,3)=c(i,3)-c(j,3) 497 498 enddo 499 500 end if 501 enddo 502 503 504200 continue 505 return 506 507 END 508c 509 subroutine smd_vlist_set(np,nl,vcutsq,c,xp,xl,point, 510 > list,cl,ct,result) 511c 512c constructs verlet pairt list 513c point(i) is an index to list() array 514c that contains all the pairs of atom i 515c In other words all the atoms paired with atom i 516c are contained in list(point(i)),...,list(point(i+1)-1) 517c np [in] number of atoms (which is also size of pointer array) 518c nl [inout] size of verlet list 519c c [in] coordinates 520c xp [in] excluded pointer 521c xl [in] excluded list 522c point [out] verlet pointer 523c list [out] verlet list 524c cl [out] list of vectors (ri-rj) 525c result [out] status of subroutine 526 implicit none 527#include "errquit.fh" 528 integer np 529 integer nl 530 double precision vcutsq 531 double precision c(np,3) 532 integer xp(np) 533 integer xl(*) 534 integer point(np) 535 integer list(nl) 536 double precision cl(nl,3) 537 double precision ct(np,3) 538 logical result 539c 540 integer i,j,k 541 integer nlist 542 integer eatm 543 double precision rij(3),rijsq,cc(1,3) 544 545 character*30 pname 546 547 pname = "smd_vlist_set" 548 549 result = .false. 550 nlist=0 551 552c write(*,*) "vcutsq=",vcutsq 553 do i=1,np-1 554 555 point(i)=nlist+1 556 if(xp(i).ne.xp(i+1))eatm=xp(i) 557 558 k = 0 559 do j=i+1,np 560 561 k = k +1 562 ct(k,1)=c(i,1)-c(j,1) 563 ct(k,2)=c(i,2)-c(j,2) 564 ct(k,3)=c(i,3)-c(j,3) 565 566 end do 567 568 569 call smd_lat_rebox(np,ct) 570 571 572 k = 0 573 do j=i+1,np 574 575 k = k + 1 576 if((xp(i).ne.xp(i+1)).and.(xl(eatm).eq.j))then 577 578 eatm=min(eatm+1,(xp(i+1)-1)) 579 580 else 581 582 rij(1)=ct(k,1) 583 rij(2)=ct(k,2) 584 rij(3)=ct(k,3) 585 586 rijsq=rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3) 587 588 589 if(rijsq.lt.vcutsq)then 590 591 nlist=nlist+1 592 593 if(nlist.gt.nl)then 594 result = .false. 595 goto 200 596 endif 597 598 list(nlist)=j 599 cl(nlist,1)=rij(1) 600 cl(nlist,2)=rij(2) 601 cl(nlist,3)=rij(3) 602 603 endif 604 605 endif 606 607 enddo 608 609 enddo 610 611 point(np)=nlist+1 612 613 nl = nlist 614 615 result = .true. 616 617200 continue 618 return 619 620 END 621 622 subroutine smd_vlist_set1(np,nl,vcutsq,c,point, 623 > list,cl,ct,result) 624c 625c constructs verlet pairt list 626c point(i) is an index to list() array 627c that contains all the pairs of atom i 628c In other words all the atoms paired with atom i 629c are contained in list(point(i)),...,list(point(i+1)-1) 630c np [in] number of atoms (which is also size of pointer array) 631c nl [inout] size of verlet list 632c c [in] coordinates 633c xp [in] excluded pointer 634c xl [in] excluded list 635c point [out] verlet pointer 636c list [out] verlet list 637c cl [out] list of vectors (ri-rj) 638c result [out] status of subroutine 639 implicit none 640#include "errquit.fh" 641 integer np 642 integer nl 643 double precision vcutsq 644 double precision c(np,3) 645 integer point(np) 646 integer list(nl) 647 double precision cl(nl,3) 648 double precision ct(np,3) 649 logical result 650c 651 integer i,j,k 652 integer nlist 653 integer eatm 654 double precision rij(3),rijsq,cc(1,3) 655 656 character*30 pname 657 658 pname = "smd_vlist_set" 659 660 result = .false. 661 nlist=0 662 663c write(*,*) "vcutsq=",vcutsq 664 do i=1,np-1 665 666 point(i)=nlist+1 667 668 k = 0 669 do j=i+1,np 670 671 k = k +1 672 ct(k,1)=c(i,1)-c(j,1) 673 ct(k,2)=c(i,2)-c(j,2) 674 ct(k,3)=c(i,3)-c(j,3) 675 676 end do 677 678 679 call smd_lat_rebox(np,ct) 680 681 682 k = 0 683 do j=i+1,np 684 685 k = k + 1 686 687 rij(1)=ct(k,1) 688 rij(2)=ct(k,2) 689 rij(3)=ct(k,3) 690 691 rijsq=rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3) 692 693 694 if(rijsq.lt.vcutsq)then 695 696 nlist=nlist+1 697 698 if(nlist.gt.nl)then 699 result = .false. 700 goto 200 701 endif 702 703 list(nlist)=j 704 cl(nlist,1)=rij(1) 705 cl(nlist,2)=rij(2) 706 cl(nlist,3)=rij(3) 707 708 endif 709 710 enddo 711 712 enddo 713 714 point(np)=nlist+1 715 716 nl = nlist 717 718 result = .true. 719 720200 continue 721 return 722 723 END 724 725 SUBROUTINE smd_vlist_test(lupdate) 726 implicit none 727#include "smd_system.fh" 728#include "mafdecls.fh" 729#include "errquit.fh" 730c 731 logical lupdate 732c 733 integer na 734 integer i_v,i_vd 735 double precision t 736 character*72 sp_vel 737 character*72 sp_vlist 738 logical result 739 character*72 tag 740 character*30 pname 741 double precision vcut,rcut 742c 743 pname = "smd_vlist_test" 744c 745c get velocity aray 746c ----------------- 747 tag = "velocity" 748 call smd_system_get_component(sp_vel,tag,result) 749 if(.not.result) 750 > call errquit( 751 > pname//'no component '//tag,0,0) 752 753 tag = "vel" 754 call smd_data_get_index(sp_vel,tag,i_v,result) 755 if(.not. result) 756 > call errquit( 757 > pname//'error getting index for'//tag,0, 0) 758 call smd_data_get_size(sp_vel,tag,na,result) 759 if(.not. result) 760 > call errquit( 761 > pname//'error getting index for'//tag,0, 0) 762 na = na/3 763c 764c get time step 765c ------------ 766 if(.not.smd_system_tstep(t)) 767 > call errquit( 768 > pname//'no time step ',0,0) 769 770c 771c get verlet displacement array 772c ----------------------------- 773 tag = "verlet_list" 774 call smd_system_get_component(sp_vlist,tag,result) 775 if(.not.result) 776 > call errquit( 777 > pname//'no component '//tag,0,0) 778 779 tag = "vlist:displacement" 780 call smd_data_get_index(sp_vlist,tag,i_vd,result) 781 if(.not. result) 782 > call errquit( 783 > pname//'error getting index for'//tag,0, 0) 784 785 call smd_cutoff_get_rcut(rcut) 786 call smd_cutoff_get_rcut_verlet(vcut) 787 788 call smd_vlist_test0(na,t,rcut,vcut, 789 > dbl_mb(i_vd),dbl_mb(i_v),lupdate) 790 791 return 792 793 END 794 795 SUBROUTINE smd_vlist_test0(natms,tstep,rcut,vcut,ivv,vvv,lupdate) 796 797 implicit none 798c 799 integer natms 800 double precision rcut,vcut 801 double precision tstep 802 logical lupdate 803 double precision ivv(natms,3) 804 double precision vvv(natms,3) 805c 806 integer i,exceed 807 808 double precision tstepsq 809 double precision dispmax,dispxsq,dispysq,dispzsq,disprsq 810 811 logical lnew 812 813 data lnew/.true./ 814 815 save lnew 816 817 tstepsq=tstep**2 818 819 if(lnew)then 820 821 lupdate=.true. 822 lnew=.false. 823 824 do i=1,natms 825 826 ivv(i,1)=0.0 827 ivv(i,2)=0.0 828 ivv(i,3)=0.0 829 830 enddo 831 832 else 833 834 lupdate=.false. 835 836 dispmax=((vcut-rcut)/2.0)**2 837 838 do i=1,natms 839 840 ivv(i,1)=ivv(i,1)+vvv(i,1) 841 ivv(i,2)=ivv(i,2)+vvv(i,2) 842 ivv(i,3)=ivv(i,3)+vvv(i,3) 843 844 enddo 845 846 exceed=0 847 848 do i=1,natms 849 850 dispxsq=ivv(i,1)**2 851 dispysq=ivv(i,2)**2 852 dispzsq=ivv(i,3)**2 853 disprsq=tstepsq*(dispxsq+dispysq+dispzsq) 854 if(disprsq.gt.dispmax) then 855 exceed=exceed+1 856 write(*,*) "verlet update disp",disprsq,dispmax 857 end if 858 if(exceed.ge.2)lupdate=.true. 859 860 enddo 861 862 if(lupdate)then 863 864 do i=1,natms 865 866 ivv(i,1)=0 867 ivv(i,2)=0 868 ivv(i,3)=0 869 870 enddo 871 872 endif 873 874 endif 875 876 return 877 878 END 879c $Id$ 880