1      subroutine argos_cafe_shks(nbonds,indexl,msb,msp,ibnd,bnd,rbnd,
2     + natoms,matoms,igan,isat,isgm,isdt,ishp,xs,ys,wgt,lseq)
3c
4      implicit none
5c
6#include "argos_cafe_common.fh"
7c
8      integer msb,msp
9      integer ibnd(msb,3)
10      real*8 bnd(msb,msp,6),rbnd(msb,2)
11      integer nbonds
12      integer indexl(nbonds)
13      integer matoms,natoms
14      integer igan(matoms),isat(matoms),isdt(matoms),ishp(matoms)
15      integer isgm(matoms),lseq(mseq)
16      real*8 xs(matoms,3),ys(matoms,3)
17      real*8 wgt(mat,mset)
18c
19      real*8 wsai,wsaj
20      logical ready,dtl
21      integer i,j,k,isb,isa,jsa,isbmax,ndtl,icset
22      integer isaglo,jsaglo
23      real*8 dmax,bs,ts,w,rsx1,rsx2,rsx3,rsi1,dif1,dif2,dif3,rsi2
24      real*8 rsi6,ff1,ff2,ff3
25c
26      character*15 filerr
27c
28#include "bitops.fh"
29c
30      niters=0
31c
32c     main iterative loop
33c
34    1 continue
35      ready=.true.
36      dtl=.false.
37      ndtl=0
38      dmax=zero
39      isbmax=0
40c
41      do 2 i=1,nbonds
42c
43      isb=indexl(i)
44c
45      isa=0
46      jsa=0
47      do 3 j=1,natoms
48      if(ibnd(isb,1).eq.igan(j)) isa=j
49      if(ibnd(isb,2).eq.igan(j)) jsa=j
50    3 continue
51c
52      if(nfhop.eq.0) then
53      icset=iset
54      else
55      icset=lseq(isgm(isa))
56      endif
57c
58      if(iand(ibnd(isb,3),icnstr).ne.0.and.
59     + bnd(isb,2,icset).gt.zero) then
60c
61      rsx1=xs(isa,1)-xs(jsa,1)
62      rsx2=xs(isa,2)-xs(jsa,2)
63      rsx3=xs(isa,3)-xs(jsa,3)
64      rbnd(isb,1)=sqrt(rsx1*rsx1+rsx2*rsx2+rsx3*rsx3)
65      rbnd(isb,2)=zero
66c      write(*,'(5i5,f12.6)') isb,isa,jsa,ishp(isa),ishp(jsa),rbnd(isb,1)
67c
68      if(iand(isdt(isa),mfixed).ne.lfixed.or.
69     + iand(isdt(jsa),mfixed).ne.lfixed) then
70c
71      bs=bnd(isb,1,icset)*bnd(isb,1,icset)
72      ts=bs*tolsha
73      wsai=one/wgt(isat(isa),icset)
74      wsaj=one/wgt(isat(jsa),icset)
75      w=wsai+wsaj
76      rsi1=bs-rsx1**2-rsx2**2-rsx3**2
77      if(abs(rsi1).gt.dmax) then
78      dmax=abs(rsi1)
79      isbmax=isb
80      endif
81      if(abs(rsi1).gt.ts) then
82      ready=.false.
83      dif1=ys(isa,1)-ys(jsa,1)
84      dif2=ys(isa,2)-ys(jsa,2)
85      dif3=ys(isa,3)-ys(jsa,3)
86      rsi2=dif1*rsx1+dif2*rsx2+dif3*rsx3
87      rsx1=dif1
88      rsx2=dif2
89      rsx3=dif3
90      if(rsi2.lt.small) then
91      rsi2=small
92      dtl=.true.
93      ndtl=isb
94      endif
95      rsi6=half*rsi1/(rsi2*w)
96      ff1=rsi6*rsx1
97      ff2=rsi6*rsx2
98      ff3=rsi6*rsx3
99      if(iand(ishp(jsa),1).ne.1) then
100      if(iand(isdt(isa),mdynam).eq.ldynam.or.
101     + iand(isdt(isa),mrestr).eq.lrestr) then
102      xs(isa,1)=xs(isa,1)+ff1*wsai
103      xs(isa,2)=xs(isa,2)+ff2*wsai
104      xs(isa,3)=xs(isa,3)+ff3*wsai
105      endif
106      endif
107      if(iand(ishp(isa),1).ne.1) then
108      if(iand(isdt(jsa),mdynam).eq.ldynam.or.
109     + iand(isdt(jsa),mrestr).eq.lrestr) then
110      xs(jsa,1)=xs(jsa,1)-ff1*wsaj
111      xs(jsa,2)=xs(jsa,2)-ff2*wsaj
112      xs(jsa,3)=xs(jsa,3)-ff3*wsaj
113      endif
114      endif
115      if(ithint) then
116      deriv(19,1)=deriv(19,1)+tstepi*tstepi*
117     + (rsx1*ff1+rsx2*ff2+rsx3*ff3)*bnd(isb,1,4)/
118     + sqrt(rsx1*rsx1+rsx2*rsx2+rsx3*rsx3)
119      endif
120      endif
121      endif
122c
123c     place dummy hydrogens (acceptor sites) onto bound heavy atoms
124c
125      if(iand(ishp(isa),1).eq.1.and.iand(ishp(jsa),1).ne.1) then
126      xs(isa,1)=xs(jsa,1)
127      xs(isa,2)=xs(jsa,2)
128      xs(isa,3)=xs(jsa,3)
129      endif
130      if(iand(ishp(isa),1).ne.1.and.iand(ishp(jsa),1).eq.1) then
131      xs(jsa,1)=xs(isa,1)
132      xs(jsa,2)=xs(isa,2)
133      xs(jsa,3)=xs(isa,3)
134      endif
135      endif
136c
137    2 continue
138c
139      niters=niters+1
140c
141      if(niters.gt.mshits) call md_abort('Too many iterations',0)
142c
143      if(dtl) then
144      isaglo=ibnd(ndtl,1)
145      jsaglo=ibnd(ndtl,2)
146      write(filerr,1000) me
147 1000 format('shake_',i3.3,'.error')
148      open(unit=16,file=filerr)
149      do 4 i=1,nbonds
150      isb=indexl(i)
151      if(iand(ibnd(isb,3),icnstr).ne.0) then
152      isa=0
153      jsa=0
154      do 5 j=1,natoms
155      if(ibnd(isb,1).eq.igan(j)) isa=j
156      if(ibnd(isb,2).eq.igan(j)) jsa=j
157    5 continue
158      if(isb.eq.ndtl) then
159      write(16,1001) igan(isa),igan(jsa),
160     + (ys(isa,k),k=1,3),(ys(jsa,k),k=1,3),
161     + (xs(isa,k),k=1,3),(xs(jsa,k),k=1,3),' < '
162 1001 format(2i7,12f8.3,a)
163      else
164      write(16,1001) igan(isa),igan(jsa),
165     + (ys(isa,k),k=1,3),(ys(jsa,k),k=1,3),
166     + (xs(isa,k),k=1,3),(xs(jsa,k),k=1,3)
167      endif
168      endif
169    4 continue
170      close(unit=16)
171      call md_abort('Deviation too large solute',isaglo)
172      endif
173c
174      if(.not.ready) goto 1
175c
176      return
177      end
178
179c $Id$
180