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