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