1      subroutine argos_cafe_polar(lself,lpbcs,
2     + xw,xwm,fw,pw,pwp,zwi,rtos,iwdt,iwz,
3     + iwfr,iwto,jwfr,jwto,xs,xsm,fs,ps,psp,zs,
4     + isga,isat,isdt,isgr,ismf,isml,isss,isq1,isq2,isq3,isgm,ishop,isz,
5     + isfr,isto,jsfr,jsto,lpbc,lstptr,lupden,eww,esw,ess,fss,esa,lseq)
6c
7c     in     log : lself           : true for box self interactions
8c     in     r*8 : xw(mwm,3,mwa)   : solvent coordinates
9c     in     r*8 : xwm(mwm,3)      : solvent molecule center of mass coordinates
10c     in/out r*8 : fw(nwm,3,mwa,2) : solvent forces
11c     out    r*8 : rtos(mwm)       : solvent distances to closes solute atom
12c     in     int : iwdt(mwm)       : solvent dynamics type
13c     out    int : iwz(mwm)        : solvent boundary type
14c     in     int : iwfr,iwto       : first and last solvent molecule i
15c     in     int : jwfr,jwto       : first and last solvent molecule j
16c     in     r*8 : xs(msa,3)       : solute atom coordinates
17c     in     r*8 : xsm(msm,3)      : solute molecule center of mass coordinates
18c     in/out r*8 : fs(msa,3,2)     : solute forces
19c     out    r*8 : zs(msm,3)       : solute virial
20c     in     int : isga(msa)       : solute global atom number
21c     in     int : isat(msa)       : solute atom type
22c     in     int : isdt(msa)       : solute dynamics type
23c     in     int : isgr(msa)       : solute charge group
24c     in     int : ismf(msa)       : solute molecule fraction
25c     in     int : isml(msa)       : solute molecule
26c     in     int : isss(msa)       : solute separation shifted scaling type
27c     in     int : isq1(msa)       : solute charge type 1
28c     in     int : isq2(msa)       : solute charge type 2
29c     in     int : isq3(msa)       : solute charge type 3
30c     out    int : isz(msa)        : solute boundary type
31c     in     int : isfr,isto       : first and last solute atom i
32c     in     int : jsfr,jsto       : first and last solute atom j
33c     in     log : lpbc            : flag for periodic boundary conditions
34c     in/out int : lstptr          : list pointer
35c     in/out r*8 : eww(mpe,2)      : mpe=1 solvent bond energy
36c                                        2 solvent angle energy
37c                                        3 solvent torsion energy
38c                                        4 solvent improper energy
39c                                        5 solvent van der Waals energy
40c                                        6 solvent electrostatic energy
41c                                        7 solvent potential energy
42c                  eswl(msf,2)     : solute-solvent van der Waals energy
43c     in/out r*8 : eswq(msf,2)     : solute-solvent electrostatic energy
44c     in/out r*8 : essl(msf,msf,2) : solute-solute van der Waals energy
45c     in/out r*8 : essq(msf,msf,2) : solute-solute electrostatic energy
46c     in/out r*8 : essr(msf,msf,2) : solute-solute reaction field energy
47c     in/out r*8 : esp(msf,2)      : solute potential energy
48c     in/out r*8 : esb(msf,2)      : solute bond energy
49c     in/out r*8 : esh(msf,2)      : solute angle energy
50c     in/out r*8 : esd(msf,2)      : solute torsion energy
51c     in/out r*8 : eso(msf,2)      : solute improper energy
52c     in     log : lupden          : update bonded energy contrbutions
53c
54c     dimensions nwm,nwa and nsa need to have been given by a call to argos_cafe_initx
55c
56      implicit none
57c
58#include "argos_cafe_common.fh"
59#include "mafdecls.fh"
60c
61      real*8 xw(mwm,3,mwa),xwm(mwm,3),fw(mwm,3,mwa,2),zwi(3,3,2)
62      real*8 xs(msa,3),xsm(msm,3),fs(msa,3,2)
63      real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2)
64      real*8 ps(msa,3,2),psp(msa,3,2,2)
65      integer iwdt(mwm),iwz(mwm),isz(msa)
66      integer isga(msa),isat(msa),isdt(msa),isgr(msa),ismf(msa)
67      integer isml(msa),isss(msa),isq1(msa),isq2(msa),isq3(msa)
68      integer isgm(msa),ishop(msa)
69      integer lseq(mseq)
70      integer iwfr,iwto,jwfr,jwto,isfr,isto,jsfr,jsto
71      integer lstptr
72      logical lself,lpbc,lupden,lpbcs,lforce
73      real*8 rtos(mwm),zs(msf,3,32)
74      real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),esa(nsa)
75      real*8 fss(msf,msf,3,2)
76c
77      integer nwloc,nsloc,nwnon,nsnon,npairs
78      integer lptr
79      integer i,j
80c
81c      print*,'CAFE forces'
82c
83      do 1 j=1,2
84      do 2 i=1,3
85      zw(1,i,j)=zwi(1,i,j)
86      zw(2,i,j)=zwi(2,i,j)
87      zw(3,i,j)=zwi(3,i,j)
88    2 continue
89    1 continue
90c
91      if(lself) then
92      jwfr=iwfr
93      jwto=iwto
94      jsfr=isfr
95      jsto=isto
96      endif
97c
98      nwloc=iwto-iwfr+1
99      if(iwfr.eq.0.or.iwto.lt.iwfr) nwloc=0
100      nwnon=jwto-jwfr+1
101      if(jwfr.eq.0.or.jwto.lt.jwfr) nwnon=0
102      nsloc=isto-isfr+1
103      if(isfr.eq.0.or.isto.lt.isfr) nsloc=0
104      nsnon=jsto-jsfr+1
105      if(jsfr.eq.0.or.jsto.lt.jsfr) nsnon=0
106c
107c     list [ list(1)=int_mb(lptr) ] with offset ndxp
108c
109c     list( 1) : index to pairlist for solvent-solvent
110c     list( 2) : index to pairlist for solute-solvent
111c     list( 3) : index to pairlist for solvent-solute
112c     list( 4) : index to pairlist for solute-solute
113c     list( 5) : index to pairlist for solute bonds
114c     list( 6) : index to pairlist for solute angles
115c     list( 7) : index to pairlist for solute dihedrals
116c     list( 8) : index to pairlist for solute impropers
117c     list( 9) : index to pairlist for solute third neighbors
118c     list(10) : index to pairlist for solute excluded pairs
119c
120c
121c     pairlists
122c     ---------
123c
124      if(lpair) call argos_pairs(lself,lpbcs,xw,xwm,iwdt,iwz,
125     + iwfr,iwto,jwfr,jwto,xs,xsm,
126     + isga,isat,isdt,isgr,isgm,ismf,isml,isss,isq1,isq2,isq3,ishop,isz,
127     + isfr,isto,jsfr,jsto,lpbc,lstptr,lseq)
128c
129      if(.not.lforce) return
130c
131c     forces
132c     ------
133c
134      ndxp=lstptr
135      lptr=i_list+ndxp+24
136c
137      if(lself.and.nwloc.gt.0) then
138      call argos_cafe_fw(iwfr,iwto,xw,fw,iwdt,int_mb(i_iwa),
139     + int_mb(i_iwq),
140     + lpbc,eww,dbl_mb(i_vdw),dbl_mb(i_chg),mbt(1),numb(1),mbp(1),
141     + int_mb(i_ibnd(1)),dbl_mb(i_bnd(1)),dbl_mb(i_rbnd(1)),
142     + mht(1),numh(1),mhp(1),
143     + int_mb(i_iang(1)),dbl_mb(i_ang(1)),dbl_mb(i_rang(1)),
144     + dbl_mb(i_rub(1)),
145     + mdt(1),numd(1),mdp(1),
146     + int_mb(i_idih(1)),dbl_mb(i_dih(1)),dbl_mb(i_rdih(1)),
147     + mit(1),numi(1),mip(1),
148     + int_mb(i_iimp(1)),dbl_mb(i_imp(1)),dbl_mb(i_rimp(1)),
149     + mtt(1),numt(1),int_mb(i_itrd(1)),
150     + mxt(1),numx(1),int_mb(i_ixcl(1)))
151      endif
152c
153c     solvent-solvent forces
154c
155      npairs=int_mb(lptr)
156      if(npairs.gt.0) then
157      call argos_cafe_fpww(xw,xwm,fw,pw,pwp,iwdt,iwfr,nwloc,lpbc,eww,
158     + dbl_mb(i_vdw),dbl_mb(i_chg),int_mb(i_iwa),int_mb(i_iwq),
159     + int_mb(i_s2i1),
160     + int_mb(lptr+1),int_mb(lptr+1+2*nwloc),int_mb(lptr+1+4*nwloc),
161     + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1),
162     + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4),
163     + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),
164     + dbl_mb(i_smr5),dbl_mb(i_smr6))
165      lptr=lptr+4*nwloc+1+npairs
166      else
167      lptr=lptr+1
168      endif
169c
170c     solute-solvent forces
171c
172      npairs=int_mb(lptr)
173      if(npairs.gt.0) then
174      call argos_cafe_fpsw(xs,xsm,fs,zs,ps,psp,isga,isat,isdt,
175     + ismf,isml,isss,
176     + isq1,isfr,nsloc,
177     + xw,xwm,fw,pw,pwp,rtos,iwdt,lpbc,lpbcs,esw,esa,
178     + dbl_mb(i_vdw),dbl_mb(i_chg),int_mb(i_iwa),int_mb(i_iwq),
179     + int_mb(i_ias),int_mb(i_s2i1),
180     + int_mb(lptr+1),int_mb(lptr+1+2*nsloc),int_mb(lptr+1+4*nsloc),
181     + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1),
182     + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4),
183     + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),dbl_mb(i_s1r6),
184     + int_mb(i_s1i1),int_mb(i_s1i2),int_mb(i_s1i3),
185     + dbl_mb(i_smr5),dbl_mb(i_smr6))
186      lptr=lptr+4*nsloc+1+npairs
187      else
188      lptr=lptr+1
189      endif
190c
191c     solvent-solute forces
192c
193      npairs=int_mb(lptr)
194      if(npairs.gt.0) then
195      call argos_cafe_fpsw(xs,xsm,fs,zs,ps,psp,isga,isat,
196     + isdt,ismf,isml,isss,
197     + isq1,jsfr,nsnon,
198     + xw,xwm,fw,pw,pwp,rtos,iwdt,lpbc,lpbcs,esw,esa,
199     + dbl_mb(i_vdw),dbl_mb(i_chg),int_mb(i_iwa),int_mb(i_iwq),
200     + int_mb(i_ias),int_mb(i_s2i1),
201     + int_mb(lptr+1),int_mb(lptr+1+2*nsnon),int_mb(lptr+1+4*nsnon),
202     + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1),
203     + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4),
204     + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),dbl_mb(i_s1r6),
205     + int_mb(i_s1i1),int_mb(i_s1i2),int_mb(i_s1i3),
206     + dbl_mb(i_smr5),dbl_mb(i_smr6))
207      lptr=lptr+4*nsnon+1+npairs
208      else
209      lptr=lptr+1
210      endif
211c
212c     solute-solute forces
213c
214      npairs=int_mb(lptr)
215      if(npairs.gt.0) then
216      call argos_cafe_fpss(xs,xsm,fs,zs,ps,psp,isga,isat,
217     + isdt,ismf,isml,isss,
218     + isq2,isq3,isfr,nsloc,lpbc,lpbcs,ess,fss,esa,
219     + dbl_mb(i_vdw),dbl_mb(i_chg),
220     + int_mb(i_ias),int_mb(i_s2i1),
221     + int_mb(lptr+1),int_mb(lptr+1+2*nsloc),int_mb(lptr+1+4*nsloc),
222     + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1),
223     + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4),
224     + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),dbl_mb(i_s1r6),
225     + int_mb(i_s1i1),int_mb(i_s1i2),int_mb(i_s1i3),int_mb(i_s1i4),
226     + int_mb(i_s1i5),dbl_mb(i_smr5),dbl_mb(i_smr6),
227     + dbl_mb(i_smr5),dbl_mb(i_smr6))
228      lptr=lptr+4*nsloc+1+npairs
229      else
230      lptr=lptr+1
231      endif
232c
233c     solute bond forces
234c
235      npairs=int_mb(lptr)
236      if(npairs.gt.0) then
237      call argos_cafe_fsb(npairs,int_mb(lptr+1),mbt(2),mbp(2),
238     + int_mb(i_ibnd(2)),dbl_mb(i_bnd(2)),dbl_mb(i_rbnd(2)),
239     + max(isto,jsto),msa,isga,isgm,ismf,isdt,isq1,dbl_mb(i_chg),
240     + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq))
241      lptr=lptr+1+npairs
242      else
243      lptr=lptr+1
244      endif
245c
246c     solute angle forces
247c
248      npairs=int_mb(lptr)
249      if(npairs.gt.0) then
250      call argos_cafe_fsh(npairs,int_mb(lptr+1),mht(2),mhp(2),
251     + int_mb(i_iang(2)),dbl_mb(i_ang(2)),dbl_mb(i_rang(2)),
252     + dbl_mb(i_rub(2)),
253     + max(isto,jsto),msa,isga,isgm,ismf,isdt,isq1,dbl_mb(i_chg),
254     + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq))
255      lptr=lptr+1+npairs
256      else
257      lptr=lptr+1
258      endif
259c
260c     solute torsion forces
261c
262      npairs=int_mb(lptr)
263      if(npairs.gt.0) then
264      call argos_cafe_fsd(npairs,int_mb(lptr+1),mdt(2),mdp(2),
265     + int_mb(i_idih(2)),dbl_mb(i_dih(2)),dbl_mb(i_rdih(2)),
266     + max(isto,jsto),msa,isga,isgm,ismf,isdt,
267     + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq))
268      lptr=lptr+1+npairs
269      else
270      lptr=lptr+1
271      endif
272c
273c     solute improper dihedral forces
274c
275      npairs=int_mb(lptr)
276      if(npairs.gt.0) then
277      call argos_cafe_fso(npairs,int_mb(lptr+1),mit(2),mip(2),
278     + int_mb(i_iimp(2)),dbl_mb(i_imp(2)),dbl_mb(i_rimp(2)),
279     + max(isto,jsto),msa,isga,isgm,ismf,isdt,
280     + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq))
281      lptr=lptr+1+npairs
282      else
283      lptr=lptr+1
284      endif
285c
286c     solute third neighbor forces
287c
288      npairs=int_mb(lptr)
289      if(npairs.gt.0) then
290      call argos_cafe_fst(npairs,int_mb(lptr+1),mtt(2),
291     + int_mb(i_itrd(2)),
292     + dbl_mb(i_vdw),dbl_mb(i_chg),
293     + max(isto,jsto),msa,isat,isga,isgm,ismf,isdt,isq3,isss,
294     + xs,fs,ess,lpbc,lpbcs,dbl_mb(i_sti),esa,int_mb(i_lseq))
295      lptr=lptr+1+npairs
296      else
297      lptr=lptr+1
298      endif
299c
300      ndxp=lptr-i_list
301c
302      do 3 j=1,2
303      do 4 i=1,3
304      zwi(1,i,j)=zw(1,i,j)
305      zwi(2,i,j)=zw(2,i,j)
306      zwi(3,i,j)=zw(3,i,j)
307    4 continue
308    3 continue
309c
310      return
311      end
312c $Id$
313