1      subroutine argos_cafe_lsw(lpbc,lpbcs,xs,isdt,isgr,isfr,isto,
2     + xwm,iwdt,iwz,
3     + iwfr,iwto,ndim,mpairs,npairs,lswjpt,lswin,lswj,list,rwx,rw)
4c
5      implicit none
6c
7#include "argos_cafe_common.fh"
8#include "bitops.fh"
9c
10      real*8 rwx(mscr,3),rw(mscr)
11      real*8 xwm(mwm,3),xs(msa,3)
12      integer ndim
13      integer lswj(*),lswjpt(ndim,2),lswin(ndim,2)
14      integer list(mwm,2),iwdt(mwm),isdt(msa),isgr(msa)
15      integer iwz(mwm)
16c
17      integer iwmfr,iwfr,iwmto,iwto,isafr,isfr,isato,isto
18      integer mpairs,npairs,isa,iwm,nlist,ilist,isa0,iisa
19c
20      integer nswi1,nswi2
21      logical lpbc,lpbcs
22c
23c     this subroutine evaluates the solute-solvent pairlist
24c
25      iwmfr=iwfr
26      iwmto=iwto
27      isafr=isfr
28      isato=isto
29c
30c      maxlst=2*mavail
31c
32      npairs=0
33c
34c     for lstype=0 the pairlist is atom based
35c
36      if(lstype.eq.0) then
37      do 1 isa=isafr,isato
38      do 3 iwm=iwmfr,iwmto
39      rwx(iwm,1)=xs(isa,1)-xwm(iwm,1)
40      rwx(iwm,2)=xs(isa,2)-xwm(iwm,2)
41      rwx(iwm,3)=xs(isa,3)-xwm(iwm,3)
42    3 continue
43      if(lpbc.or.lpbcs)
44     + call argos_cafe_pbc(1,rwx,mscr,rwx,mscr,0,iwmfr,iwmto)
45      do 8 iwm=iwmfr,iwmto
46      rw(iwm)=rwx(iwm,1)**2+rwx(iwm,2)**2+rwx(iwm,3)**2
47    8 continue
48      do 4 iwm=iwmfr,iwmto
49      list(iwm,1)=0
50      list(iwm,2)=0
51      if(rw(iwm).lt.rshrt2.and.
52     + (iand(iwdt(iwm),mdynam).eq.ldynam.or.
53     + iand(isdt(isa),mdynam).eq.ldynam)) list(iwm,1)=1
54      if(rw(iwm).lt.rrest2.and.
55     + (iand(iwdt(iwm),mrestr).eq.lrestr.or.
56     + iand(isdt(isa),mrestr).eq.lrestr) ) list(iwm,1)=1
57    4 continue
58      if(npsw.eq.2) then
59      do 401 iwm=iwmfr,iwmto
60      if(rw(iwm).lt.rlong2.and.
61     + (iand(iwdt(iwm),mdynam).eq.ldynam.or.
62     + iand(isdt(isa),mdynam).eq.ldynam)) list(iwm,2)=1
63      if(list(iwm,1).eq.1) list(iwm,2)=0
64  401 continue
65      endif
66c
67c     short range pairlist
68c
69      nlist=0
70      do 402 iwm=iwmfr,iwmto
71      if(list(iwm,1).eq.1) then
72      nlist=nlist+1
73      list(nlist,1)=iwm
74      endif
75  402 continue
76      if(npairs+nlist.gt.mpairs)
77     + call md_abort('Insufficient memory for pairlist',0)
78      do 5 ilist=1,nlist
79      lswj(npairs+ilist)=list(ilist,1)
80    5 continue
81      lswjpt(isa-isafr+1,1)=npairs+1
82      npairs=npairs+nlist
83      lssw=lssw+nlist
84      lswin(isa-isafr+1,1)=nlist
85c
86c     long range pairlist
87c
88      if(npsw.eq.2) then
89      nlist=0
90      do 8402 iwm=iwmfr,iwmto
91      if(list(iwm,2).eq.1) then
92      nlist=nlist+1
93      list(nlist,2)=iwm
94      endif
95 8402 continue
96      if(npairs+nlist.gt.mpairs)
97     + call md_abort('Insufficient memory for pairlist',0)
98      do 8005 ilist=1,nlist
99      lswj(npairs+ilist)=list(ilist,2)
100 8005 continue
101      lswjpt(isa-isafr+1,2)=npairs+1
102      npairs=npairs+nlist
103      llsw=llsw+nlist
104      lswin(isa-isafr+1,2)=nlist
105      endif
106    1 continue
107      endif
108c
109c     for lstype=1 the pairlist is charge group based
110c         lstype=1 the pairlist is segment based
111c
112      if(lstype.eq.1.or.lstype.eq.2) then
113      isa=isafr
114      isa0=isafr
115   21 continue
116      nlist=0
117      do 15 iwm=iwmfr,iwmto
118      list(iwm,1)=0
119      list(iwm,2)=0
120   15 continue
121   10 continue
122      do 11 iwm=iwmfr,iwmto
123      rwx(iwm,1)=xs(isa,1)-xwm(iwm,1)
124      rwx(iwm,2)=xs(isa,2)-xwm(iwm,2)
125      rwx(iwm,3)=xs(isa,3)-xwm(iwm,3)
126   11 continue
127      if(lpbc.or.lpbcs)
128     + call argos_cafe_pbc(1,rwx,mscr,rwx,mscr,0,iwmfr,iwmto)
129      do 13 iwm=iwmfr,iwmto
130      rw(iwm)=rwx(iwm,1)**2+rwx(iwm,2)**2+rwx(iwm,3)**2
131   13 continue
132      do 9012 iwm=iwmfr,iwmto
133      if(iand(iwdt(iwm),mdynam).eq.ldynam.and.
134     + rw(iwm).lt.rshrt2) list(iwm,1)=1
135      if(iand(iwdt(iwm),mrestr).eq.lrestr.and.
136     + rw(iwm).lt.rrest2) list(iwm,1)=1
137 9012 continue
138      if(iand(isdt(isa),mdynam).eq.ldynam) then
139      do 9013 iwm=iwmfr,iwmto
140      if(rw(iwm).lt.rshrt2) list(iwm,1)=1
141 9013 continue
142      endif
143      if(iand(isdt(isa),mrestr).eq.lrestr) then
144      do 9014 iwm=iwmfr,iwmto
145      if(rw(iwm).lt.rrest2) list(iwm,1)=1
146 9014 continue
147      endif
148      if(npsw.eq.2) then
149      do 9812 iwm=iwmfr,iwmto
150      if(iand(iwdt(iwm),mdynam).eq.ldynam.and.
151     + rw(iwm).lt.rlong2) list(iwm,2)=1
152 9812 continue
153      if(iand(isdt(isa),mdynam).eq.ldynam) then
154      do 9813 iwm=iwmfr,iwmto
155      if(rw(iwm).lt.rlong2) list(iwm,2)=1
156 9813 continue
157      endif
158      endif
159      do 9815 iwm=iwmfr,iwmto
160      if(list(iwm,1).eq.1) list(iwm,2)=0
161 9815 continue
162      if(isa.eq.isato) goto 16
163      if(lstype.eq.1) then
164      if(isgr(isa+1).eq.isgr(isa)) then
165      isa=isa+1
166      goto 10
167      endif
168      endif
169c      if(lstype.eq.2) then
170c      if(isl(isa+1,lssgm).eq.isl(isa,lssgm)) then
171c      isa=isa+1
172c      goto 10
173c      endif
174c      endif
175   16 continue
176c
177c     at this point the partial list for a solute charge group or
178c     segment with all solvent has been evaluated
179c
180c     short range pairlist
181c
182      nlist=0
183      do 17 iwm=iwmfr,iwmto
184      if(list(iwm,1).eq.1) then
185      nlist=nlist+1
186      list(nlist,1)=iwm
187      endif
188   17 continue
189      do 18 iisa=isa0,isa
190      lswjpt(iisa-isafr+1,1)=npairs+1
191      lswin(iisa-isafr+1,1)=nlist
192   18 continue
193      if(npairs+nlist.gt.mpairs)
194     + call md_abort('Insufficient memory for pairlist',0)
195      do 19 ilist=1,nlist
196      lswj(npairs+ilist)=list(ilist,1)
197   19 continue
198      npairs=npairs+nlist
199      lssw=lssw+nlist
200c
201c     long range pairlist
202c
203      if(npsw.eq.2) then
204      nlist=0
205      do 8017 iwm=iwmfr,iwmto
206      if(list(iwm,2).eq.1) then
207      nlist=nlist+1
208      list(nlist,2)=iwm
209      endif
210 8017 continue
211      do 8018 iisa=isa0,isa
212      lswjpt(iisa-isafr+1,2)=npairs+1
213      lswin(iisa-isafr+1,2)=nlist
214 8018 continue
215      if(npairs+nlist.gt.mpairs)
216     + call md_abort('Insufficient memory for pairlist',0)
217      do 8019 ilist=1,nlist
218      lswj(npairs+ilist)=list(ilist,2)
219 8019 continue
220      npairs=npairs+nlist
221      llsw=llsw+nlist
222      endif
223c
224      isa=isa+1
225      isa0=isa
226      if(isa.le.isato) goto 21
227      endif
228      nswi1=0
229      nswi2=0
230      do 9191 isa=isafr,isato
231      nswi1=nswi1+lswin(isa-isafr+1,1)
232      if(npsw.eq.2) nswi2=nswi2+lswin(isa-isafr+1,2)
233 9191 continue
234c
235      return
236      end
237c $Id$
238