1      subroutine argos_cafe_shkw(nwmloc,iwgm,iwdt,xw,yw,mwb,
2     + nwb,nbp,ibnd,bnd,iwat,wgt,rwx,rwi1,rwi2,rwi6,f,swt)
3c
4      implicit none
5c
6#include "argos_cafe_common.fh"
7c
8      integer nwmloc,iwgm(mwm),iwdt(mwm)
9      real*8 xw(mwm,3,mwa),yw(mwm,3,mwa)
10      integer mwb,nwb,nbp
11      integer ibnd(mwb,3),iwat(mwa)
12      real*8 bnd(mwb,nbp,6),wgt(mat,mset)
13      real*8 rwx(mscr,3),rwi1(mscr),rwi2(mscr),rwi6(mscr),f(mscr)
14      real*8 swt(mwm,3,mwb)
15c
16      logical ready
17c
18      integer ix,iwb,iwm,iwa,jwa,nshake,ndtl
19      real*8 bw,tw,w,dif,tisum,cdwb4,ydif1,ydif2,ydif3
20      real*8 ydist,ep2sum,cdwb2,ep3sum,cdwb3,wiwa,wjwa
21c
22#include "bitops.fh"
23c
24      if(ithint.and.ith(7)) then
25      do 6 ix=1,3
26      do 5 iwb=1,nwb
27      do 4 iwm=1,nwmloc
28      swt(iwm,ix,iwb)=zero
29    4 continue
30    5 continue
31    6 continue
32      endif
33c
34      niterw=0
35    1 continue
36      ready=.true.
37      do 2 iwb=1,nwb
38      if(iand(ibnd(iwb,3),icnstr).eq.0) goto 2
39      iwa=ibnd(iwb,1)
40      jwa=ibnd(iwb,2)
41      wiwa=one/wgt(iwat(iwa),iset)
42      wjwa=one/wgt(iwat(jwa),iset)
43      bw=bnd(iwb,1,iset)*bnd(iwb,1,iset)
44      tw=bw*tolsha
45      w=wiwa+wjwa
46      nshake=0
47      do 3 iwm=1,nwmloc
48      rwx(iwm,1)=xw(iwm,1,iwa)-xw(iwm,1,jwa)
49      rwx(iwm,2)=xw(iwm,2,iwa)-xw(iwm,2,jwa)
50      rwx(iwm,3)=xw(iwm,3,iwa)-xw(iwm,3,jwa)
51      rwi1(iwm)=bw-rwx(iwm,1)**2-rwx(iwm,2)**2-rwx(iwm,3)**2
52      if(abs(rwi1(iwm)).gt.tw) nshake=nshake+1
53    3 continue
54      if(nshake.gt.0) then
55      ready=.false.
56      do 7 iwm=1,nwmloc
57      dif=yw(iwm,1,iwa)-yw(iwm,1,jwa)
58      rwi2(iwm)=dif*rwx(iwm,1)
59      rwx(iwm,1)=dif
60      dif=yw(iwm,2,iwa)-yw(iwm,2,jwa)
61      rwi2(iwm)=rwi2(iwm)+dif*rwx(iwm,2)
62      rwx(iwm,2)=dif
63      dif=yw(iwm,3,iwa)-yw(iwm,3,jwa)
64      rwi2(iwm)=rwi2(iwm)+dif*rwx(iwm,3)
65      rwx(iwm,3)=dif
66    7 continue
67      ndtl=0
68      do 8 iwm=1,nwmloc
69      if(rwi2(iwm).lt.small) ndtl=iwgm(iwm)
70    8 continue
71      if(ndtl.gt.0)
72     + call md_abort('Deviation too large for solvent',ndtl)
73      do 9 iwm=1,nwmloc
74      rwi6(iwm)=half*rwi1(iwm)/(rwi2(iwm)*w)
75    9 continue
76c
77      if(ithint.and.ith(7)) then
78      do 12 ix=1,3
79      do 11 iwm=1,nwmloc
80      swt(iwm,ix,iwb)=swt(iwm,ix,iwb)+rwi6(iwm)*rwx(iwm,ix)
81   11 continue
82   12 continue
83      endif
84c
85      do 10 iwm=1,nwmloc
86      f(iwm)=rwi6(iwm)*rwx(iwm,1)
87      xw(iwm,1,iwa)=xw(iwm,1,iwa)+f(iwm)*wiwa
88      xw(iwm,1,jwa)=xw(iwm,1,jwa)-f(iwm)*wjwa
89      f(iwm)=rwi6(iwm)*rwx(iwm,2)
90      xw(iwm,2,iwa)=xw(iwm,2,iwa)+f(iwm)*wiwa
91      xw(iwm,2,jwa)=xw(iwm,2,jwa)-f(iwm)*wjwa
92      f(iwm)=rwi6(iwm)*rwx(iwm,3)
93      xw(iwm,3,iwa)=xw(iwm,3,iwa)+f(iwm)*wiwa
94      xw(iwm,3,jwa)=xw(iwm,3,jwa)-f(iwm)*wjwa
95   10 continue
96      endif
97    2 continue
98      niterw=niterw+1
99      if(niterw.gt.mshitw)
100     + call md_abort('Too many shake iterations',niterw)
101      if(.not.ready) goto 1
102c
103      do 24 iwm=1,nwmloc
104      if(iand(iwdt(iwm),mfixed).eq.lfixed) then
105      do 25 iwa=1,mwa
106      xw(iwm,1,iwa)=yw(iwm,1,iwa)
107      xw(iwm,2,iwa)=yw(iwm,2,iwa)
108      xw(iwm,3,iwa)=yw(iwm,3,iwa)
109   25 continue
110      if(ithint.and.ith(7)) then
111      do 26 iwb=1,nwb
112      swt(iwm,1,iwb)=zero
113      swt(iwm,2,iwb)=zero
114      swt(iwm,3,iwb)=zero
115   26 continue
116      endif
117      endif
118   24 continue
119c
120c      if(nwr.gt.0) then
121c      if(iropt.eq.2.or.(iropt.eq.1.and.forest.gt.tiny)) then
122c      do 21 iwr=1,nwr
123c      cmx1=zero
124c      cmx2=zero
125c      cmx3=zero
126c      do 22 iwa=1,mwa
127c      cmx1=cmx1+xw(idwr(iwr),1,iwa)*wwaf(iwa)
128c      cmx2=cmx2+xw(idwr(iwr),2,iwa)*wwaf(iwa)
129c      cmx3=cmx3+xw(idwr(iwr),3,iwa)*wwaf(iwa)
130c   22 continue
131c      do 23 iwa=1,mwa
132c      xw(idwr(iwr),1,iwa)=xw(idwr(iwr),1,iwa)+xwr(iwr,1)-cmx1
133c      xw(idwr(iwr),2,iwa)=xw(idwr(iwr),2,iwa)+xwr(iwr,2)-cmx2
134c      xw(idwr(iwr),3,iwa)=xw(idwr(iwr),3,iwa)+xwr(iwr,3)-cmx3
135c   23 continue
136c   21 continue
137c      endif
138c      endif
139c
140      if(ithint.and.ith(7)) then
141      tisum=zero
142      do 14 iwb=1,nwb
143      iwa=ibnd(iwb,1)
144      jwa=ibnd(iwb,2)
145      cdwb4=bnd(iwb,1,4)
146      do 15 iwm=1,nwmloc
147      ydif1=yw(iwm,1,iwa)-yw(iwm,1,jwa)
148      ydif2=yw(iwm,2,iwa)-yw(iwm,2,jwa)
149      ydif3=yw(iwm,3,iwa)-yw(iwm,3,jwa)
150      ydist=sqrt(ydif1**2+ydif2**2+ydif3**2)
151      tisum=tisum+(swt(iwm,1,iwb)*ydif1+swt(iwm,2,iwb)*ydif2+
152     + swt(iwm,3,iwb)*ydif3)*cdwb4/ydist
153   15 continue
154   14 continue
155      deriv(7,1)=tisum*tstepi*tstepi
156      endif
157      if(ip2(5)) then
158      ep2sum=zero
159      do 16 iwb=1,nwb
160      iwa=ibnd(iwb,1)
161      jwa=ibnd(iwb,2)
162      cdwb2=bnd(iwb,1,2)-bnd(iwb,1,1)
163      do 17 iwm=1,nwmloc
164      ydif1=yw(iwm,1,iwa)-yw(iwm,1,jwa)
165      ydif2=yw(iwm,2,iwa)-yw(iwm,2,jwa)
166      ydif3=yw(iwm,3,iwa)-yw(iwm,3,jwa)
167      ydist=sqrt(ydif1**2+ydif2**2+ydif3**2)
168      ep2sum=ep2sum+(swt(iwm,1,iwb)*ydif1+swt(iwm,2,iwb)*ydif2+
169     + swt(iwm,3,iwb)*ydif3)*cdwb2/ydist
170   17 continue
171   16 continue
172      ep2(1)=ep2(1)+ep2sum*tstepi*tstepi
173      endif
174      if(ip3(5)) then
175      ep3sum=zero
176      do 18 iwb=1,nwb
177      iwa=ibnd(iwb,1)
178      jwa=ibnd(iwb,2)
179      cdwb3=bnd(iwb,1,3)-bnd(iwb,1,1)
180      do 19 iwm=1,nwmloc
181      ydif1=yw(iwm,1,iwa)-yw(iwm,1,jwa)
182      ydif2=yw(iwm,2,iwa)-yw(iwm,2,jwa)
183      ydif3=yw(iwm,3,iwa)-yw(iwm,3,jwa)
184      ydist=sqrt(ydif1**2+ydif2**2+ydif3**2)
185      ep3sum=ep3sum+(swt(iwm,1,iwb)*ydif1+swt(iwm,2,iwb)*ydif2+
186     + swt(iwm,3,iwb)*ydif3)*cdwb3/ydist
187   19 continue
188   18 continue
189      ep3(1)=ep3(1)+ep3sum*tstepi*tstepi
190      endif
191c
192      return
193      end
194c $Id$
195