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