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