1
2      subroutine smd_exlist_init_system()
3      implicit none
4#include "errquit.fh"
5#include "inp.fh"
6#include "mafdecls.fh"
7#include "rtdb.fh"
8#include "util.fh"
9#include "global.fh"
10c
11      character*32 sp_exlist,sp_atom
12      character*32 tag,pname
13      logical result
14
15      pname = "smd_exlist_init_system"
16c
17      tag = "atom"
18      call smd_system_get_component(sp_atom,tag,result)
19      if(.not.result)
20     >  call errquit(
21     >       pname//'no component '//tag,0,0)
22
23      tag = "excl_list"
24      call smd_system_get_component(sp_exlist,tag,result)
25      if(.not.result)
26     >  call errquit(
27     >       pname//'no component '//tag,0,0)
28
29
30      call smd_exlist_init(sp_exlist,sp_atom,result)
31c
32      if(.not.result) then
33       tag = "excl_list"
34       call smd_system_unset_component(tag)
35      end if
36
37      return
38      end
39
40      subroutine smd_exlist_init(sp_exlist,sp_atom,result)
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_exlist
50      character*(*) sp_atom
51      logical result
52c
53      character*32 pname
54      character*80 tag
55      integer np,nl
56      integer i_p,i_l,i_ir
57      integer h_l
58      integer i_list
59      integer i
60c
61      pname = "smd_exlist_init"
62c
63c      write(*,*) "in "//pname
64c
65c     get total number of atoms
66c     -------------------------
67      call smd_atom_ntot(sp_atom,np)
68c
69c     gestimate the size of pair list
70c     ------------------------------
71      nl =  min( np*200, ma_inquire_avail(MT_INT))
72c
73c     create pointer memory
74c     ---------------------
75      call smd_namespace_create(sp_exlist)
76      tag = "exlist:pointer"
77      call smd_data_create(sp_exlist,"exlist:pointer",np,MT_INT)
78      call smd_data_get_index(sp_exlist,tag,i_p,result)
79      if(.not. result)
80     >  call errquit(
81     >       pname//'error getting index for'//tag,0, RTDB_ERR)
82
83c
84c    create temporary scratch array for list since
85c    we do not know the size yet
86c    ---------------------------------------------
87      if(.not.ma_push_get(mt_int,nl,'tmp l',h_l,i_l))
88     + call errquit(pname//'Failed to allocate memory for tmp l',
89     + nl, MA_ERR)
90
91c
92c     get atom residue index for exclusion criteria
93c     ---------------------------------------------
94      tag = "atom:resid"
95      call smd_data_get_index(sp_atom,tag,i_ir,result)
96      if(.not. result)
97     >  call errquit(
98     >       pname//'error getting index for'//tag,0, RTDB_ERR)
99
100
101      call smd_exlist_set(nl,
102     +                    int_mb(i_l),
103     +                    np,
104     +                    int_mb(i_p),
105     +                    int_mb(i_ir))
106
107c
108c     create list memory
109c     nl now contains the actual size
110c     -------------------------------
111
112      if(nl.eq.0) then
113        result = .false.
114        call smd_namespace_destroy(sp_exlist)
115        goto 200
116      end if
117c
118c      write(*,*) "excluded size list",nl
119      tag = "exlist:list"
120      call smd_data_create(sp_exlist,tag,nl,MT_INT)
121      call smd_data_get_index(sp_exlist,tag,i_list,result)
122      if(.not. result)
123     >  call errquit(
124     >       pname//'error getting index for'//tag,0, RTDB_ERR)
125
126
127      do i=1,nl
128       int_mb(i_list+i-1) = int_mb(i_l+i-1)
129      end do
130
131      if(.not.ma_pop_stack(h_l))
132     & call errquit(pname//'Failed to deallocate stack h_l',nl,
133     &       MA_ERR)
134
135200   continue
136      return
137      end
138c
139c
140      subroutine smd_exlist_set(nl,list,natms,point,resid)
141c
142c     constructs excluded pairt list
143c     point(i) is an index to list() array
144c     that contains all the pairs of atom i
145c     In other words all the atoms paired with atom i
146c     are contained in list(point(i)),...,list(point(i+1)-1)
147c
148      implicit none
149#include "errquit.fh"
150      integer nl,natms
151      integer list(nl)
152      integer point(natms)
153      integer resid(natms)
154c
155      integer i,j
156      integer nlist
157
158      character*30 pname
159
160      pname = "smd_exlist_set"
161
162      nlist = 0
163      do i=1,natms-1
164
165       point(i)=nlist+1
166
167       do j=i+1,natms
168
169         if(resid(i).eq.resid(j))then
170
171          nlist=nlist+1
172
173          if(nlist.gt.nl)
174     +     call errquit(pname//'exceeded list dimensions',
175     +     nl, 0)
176
177          list(nlist)=j
178
179         endif
180
181       enddo
182      enddo
183c
184c     we need to set this to mark
185c     the end of pair list index for natms-1
186c     it would be used as point(natms)-1
187c     -------------------------------------
188      point(natms)=nlist+1
189
190      nl = nlist
191
192      return
193
194      END
195c $Id$
196