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