1 subroutine argos_cafe_fw(iwfr,iwto,xw,fw,iwdt,iwatm,iwq, 2 + lpbc,eww,vdw,chg, 3 + mwb,nwb,nbp,ibnd,bnd,rbnd,mwh,nwh,nhp,iang,ang,rang,rub, 4 + mwd,nwd,ndp,idih,dih,rdih,mwo,nwo,nop,iimp,dimp,rimp, 5 + mwt,nwt,idwt,mwn,nwn,idwn) 6c 7c $Id$ 8c 9 implicit none 10c 11#include "argos_cafe_common.fh" 12c 13 integer iwfr,iwto 14 integer mwb,mwh,mwd,mwo,nbp,nhp,ndp,nop,mwt,mwn 15 integer nwb,nwh,nwd,nwo,nwt,nwn 16 real*8 xw(mwm,3,mwa),fw(mwm,3,mwa,2) 17 integer iwdt(mwm),iwq(mwa),iwatm(mwa) 18 logical lpbc 19 real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset) 20 integer ibnd(mwb,3),iang(mwh,4),idih(mwd,5),iimp(mwo,5) 21 real*8 bnd(mwb,nbp,mset),ang(mwh,nhp,mset) 22 real*8 dih(mwd,ndp,mset),dimp(mwo,nop,mset) 23 real*8 rbnd(mwb,2),rang(mwh,2),rub(mwh,2),rdih(mwd,2),rimp(mwo,2) 24c 25c real*8 ca6(mat,mat,6),ca12(mat,mat,6) 26c real*8 cb6(mat,mat,6),cb12(mat,mat,6) 27c integer iwl(mwm,miw2), 28c 29 integer idwt(0:mwt,2),idwn(0:mwn,2) 30c 31c real*8 cdwb(mwb,6),ddwb(mwb,6) 32c integer iwbs(mwb),idwb(mwb),jdwb(mwb),iwatm(mwa) 33c real*8 cdwh(mwh,6),ddwh(mwh,6) 34c integer idwh(mwh),jdwh(mwh),kdwh(mwh) 35c real*8 cdwd(mwd,6),ddwd(mwd,6),edwd(mwd,6) 36c integer idwd(mwd),jdwd(mwd),kdwd(mwd),ldwd(mwd) 37c real*8 cdwo(mwo,6),ddwo(mwo,6) 38c integer idwo(mwo),jdwo(mwo),kdwo(mwo),ldwo(mwo) 39c real*8 uwb(mwb),uwh(mwh),uwd(mwd),uwo(mwo) 40c 41 integer iwb,iwa,jwa,iwm,iwh,kwa,iwd,lwa,iwo,iwt,iwn 42 real*8 bond,for,rwx1,rwx2,rwx3,rww,rwwi,dbond,dfor,dfw1,dfw2,dfw3 43 real*8 angle,xwij1,xwij2,xwij3,xwkj1,xwkj2,xwkj3,rwij2,rwij2i 44 real*8 rwkj2,rwkj2i,cphi,phi,dangle,sphi,rmul 45 real*8 xwkl1,xwkl2,xwkl3,xwik1,xwik2,xwik3,xwjl1,xwjl2,xwjl3 46 real*8 xm1,xm2,xm3,xn1,xn2,xn3,rm2i,rn2i,rmni,s,rpa 47 real*8 xd1,xd2,xd3,xe1,xe2,xe3,dfwi1,dfwi2,dfwi3 48 real*8 dfwj1,dfwj2,dfwj3,dfwk1,dfwk2,dfwk3,dfwl1,dfwl2,dfwl3 49 real*8 danglep,c6p1,c12p1,c6p2,c12p2,qip1,qjp1,qip2,qjp2 50 real*8 c6,c12,c6t,c12t,qit,qjt,cf6,cf12,qi,qj,q,qp1,qp2 51 real*8 ep2l,ep3l,ep2q,ep3q,rxx,rxy,rxz,r2,r2i,r1i,r6i,dfw 52 real*8 rwikji,sphii,qij,rwi,ferfc,fderfc,eww(mpe,2) 53 real*8 etermq,eterml,eub 54c 55#include "argos_cafe_funcs_dec.fh" 56#include "bitops_decls.fh" 57#include "argos_cafe_funcs_sfn.fh" 58#include "bitops_funcs.fh" 59c 60 c6t=zero 61 c12t=zero 62 qit=zero 63 qjt=zero 64 qp1=zero 65 qp2=zero 66c 67 do 10 iwb=1,nwb 68 if(iand(ibnd(iwb,3),icnstr).eq.0) then 69 iwa=ibnd(iwb,1) 70 jwa=ibnd(iwb,2) 71 bond=bnd(iwb,1,iset) 72 for=bnd(iwb,2,iset) 73 rbnd(iwb,2)=zero 74 do 20 iwm=iwfr,iwto 75 rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 76 rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 77 rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 78 rww=sqrt(rwx1**2+rwx2**2+rwx3**2) 79 if(rww.lt.tiny) then 80 rwwi=one 81 else 82 rwwi=one/rww 83 endif 84 dbond=rww-bond 85 if(iand(iwdt(iwm),mdynam).eq.ldynam) 86 + rbnd(iwb,2)=rbnd(iwb,2)+half*for*(rww-bond)**2 87 dfor=for*dbond*rwwi 88 dfw1=dfor*rwx1 89 dfw2=dfor*rwx2 90 dfw3=dfor*rwx3 91 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1 92 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1 93 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2 94 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2 95 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3 96 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3 97 if(ip2(6)) 98 + ep2(1)=ep2(1)+half*bnd(iwb,2,2)*(rww-bnd(iwb,1,2))**2 99 if(ip3(6)) 100 + ep3(1)=ep3(1)+half*bnd(iwb,2,3)*(rww-bnd(iwb,1,3))**2 101 if(ith(6)) then 102 deriv(6,1)=deriv(6,1)+ 103 + dbond*(half*dbond*bnd(iwb,2,4)-for*bnd(iwb,1,4)) 104 endif 105 20 continue 106 eww(1,1)=eww(1,1)+rbnd(iwb,2) 107 if(ip2(6)) ep2(1)=ep2(1)-rbnd(iwb,2) 108 if(ip3(6)) ep3(1)=ep3(1)-rbnd(iwb,2) 109 endif 110 if(ipme.ne.0) then 111 iwa=ibnd(iwb,1) 112 jwa=ibnd(iwb,2) 113 qij=chg(iwq(iwa),1,iset)*chg(iwq(jwa),1,iset) 114 do 21 iwm=iwfr,iwto 115 rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 116 rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 117 rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 118 rww=sqrt(rwx1**2+rwx2**2+rwx3**2) 119 rwi=one/rww 120 ferfc=one-erfc(ealpha*rww) 121 fderfc=-(ealpha*derfc(ealpha*rww)) 122 epmecw=epmecw-ferfc*qij*rwi 123 eww(9,1)=eww(9,1)-ferfc*qij*rwi 124 dfor=-(qij*rwi*rwi*(ferfc*rwi-fderfc)) 125 dfw1=dfor*rwx1 126 dfw2=dfor*rwx2 127 dfw3=dfor*rwx3 128 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1 129 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1 130 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2 131 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2 132 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3 133 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3 134 vpmeb(1)=vpmeb(1)+dfw1*rwx1 135 vpmeb(2)=vpmeb(2)+dfw2*rwx1 136 vpmeb(3)=vpmeb(3)+dfw3*rwx1 137 vpmeb(4)=vpmeb(4)+dfw2*rwx2 138 vpmeb(5)=vpmeb(5)+dfw3*rwx2 139 vpmeb(6)=vpmeb(6)+dfw3*rwx3 140 21 continue 141 endif 142 10 continue 143 do 40 iwh=1,nwh 144 iwa=iang(iwh,1) 145 jwa=iang(iwh,2) 146 kwa=iang(iwh,3) 147 angle=ang(iwh,1,iset) 148 for=ang(iwh,2,iset) 149 rang(iwh,2)=zero 150 do 50 iwm=iwfr,iwto 151 xwij1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 152 xwij2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 153 xwij3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 154 xwkj1=xw(iwm,1,kwa)-xw(iwm,1,jwa) 155 xwkj2=xw(iwm,2,kwa)-xw(iwm,2,jwa) 156 xwkj3=xw(iwm,3,kwa)-xw(iwm,3,jwa) 157 rwij2=xwij1**2+xwij2**2+xwij3**2 158 rwkj2=xwkj1**2+xwkj2**2+xwkj3**2 159 rwij2i=one/rwij2 160 rwkj2i=one/rwkj2 161 rwikji=one/sqrt(rwij2*rwkj2) 162 cphi=rwikji*(xwij1*xwkj1+xwij2*xwkj2+xwij3*xwkj3) 163 if(cphi.lt.-one) cphi=-one 164 if(cphi.gt. one) cphi= one 165 phi=acos(cphi) 166 dangle=phi-angle 167 if(iand(iwdt(iwm),mdynam).eq.ldynam) 168 + rang(iwh,2)=rang(iwh,2)+half*for*dangle*dangle 169 sphi=sin(phi) 170 if(abs(sphi).lt.small) sphi=small 171 dfor=for*dangle/sphi 172 dfw1=dfor*(xwkj1*rwikji-xwij1*rwij2i*cphi) 173 dfw2=dfor*(xwkj2*rwikji-xwij2*rwij2i*cphi) 174 dfw3=dfor*(xwkj3*rwikji-xwij3*rwij2i*cphi) 175 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)+dfw1 176 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw1 177 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)+dfw2 178 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw2 179 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)+dfw3 180 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw3 181 dfw1=dfor*(xwij1*rwikji-xwkj1*rwkj2i*cphi) 182 dfw2=dfor*(xwij2*rwikji-xwkj2*rwkj2i*cphi) 183 dfw3=dfor*(xwij3*rwikji-xwkj3*rwkj2i*cphi) 184 fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)+dfw1 185 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw1 186 fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)+dfw2 187 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw2 188 fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)+dfw3 189 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw3 190 if(ip2(8)) 191 + ep2(1)=ep2(1)+half*ang(iwh,2,2)*(phi-ang(iwh,1,2))**2 192 if(ip3(8)) 193 + ep3(1)=ep3(1)+half*ang(iwh,2,3)*(phi-ang(iwh,1,3))**2 194 if(ith(8)) then 195 deriv(8,1)=deriv(8,1)+ 196 + dangle*(half*dangle*ang(iwh,2,4)-for*ang(iwh,1,4)) 197 endif 198 50 continue 199 eww(2,1)=eww(2,1)+rang(iwh,2) 200 if(ip2(8)) ep2(1)=ep2(1)-rang(iwh,2) 201 if(ip3(8)) ep3(1)=ep3(1)-rang(iwh,2) 202 if(ipme.ne.0) then 203 iwa=iang(iwh,1) 204 jwa=iang(iwh,3) 205 qij=chg(iwq(iwa),1,iset)*chg(iwq(jwa),1,iset) 206 do 41 iwm=iwfr,iwto 207 rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 208 rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 209 rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 210 rww=sqrt(rwx1**2+rwx2**2+rwx3**2) 211 rwi=one/rww 212 ferfc=one-erfc(ealpha*rww) 213 fderfc=-(ealpha*derfc(ealpha*rww)) 214 epmecw=epmecw-ferfc*qij*rwi 215 eww(9,1)=eww(9,1)-ferfc*qij*rwi 216 dfor=-(qij*rwi*rwi*(ferfc*rwi-fderfc)) 217 dfw1=dfor*rwx1 218 dfw2=dfor*rwx2 219 dfw3=dfor*rwx3 220 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1 221 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1 222 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2 223 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2 224 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3 225 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3 226 vpmeb(1)=vpmeb(1)+dfw1*rwx1 227 vpmeb(2)=vpmeb(2)+dfw2*rwx1 228 vpmeb(3)=vpmeb(3)+dfw3*rwx1 229 vpmeb(4)=vpmeb(4)+dfw2*rwx2 230 vpmeb(5)=vpmeb(5)+dfw3*rwx2 231 vpmeb(6)=vpmeb(6)+dfw3*rwx3 232 41 continue 233 endif 234 40 continue 235 if(iffld.eq.2) then 236 do 1140 iwh=1,nwh 237 iwa=iang(iwh,1) 238 kwa=iang(iwh,3) 239 bond=ang(iwb,3,iset) 240 for=ang(iwb,4,iset) 241 eub=zero 242 do 150 iwm=iwfr,iwto 243 rwx1=xw(iwm,1,iwa)-xw(iwm,1,kwa) 244 rwx2=xw(iwm,2,iwa)-xw(iwm,2,kwa) 245 rwx3=xw(iwm,3,iwa)-xw(iwm,3,kwa) 246 rww=sqrt(rwx1**2+rwx2**2+rwx3**2) 247 if(rww.lt.tiny) then 248 rwwi=one 249 else 250 rwwi=one/rww 251 endif 252 dbond=rww-bond 253 if(iand(iwdt(iwm),mdynam).eq.ldynam) 254 + eub=eub+half*for*(rww-bond)**2 255 dfor=for*dbond*rwwi 256 dfw1=dfor*rwx1 257 dfw2=dfor*rwx2 258 dfw3=dfor*rwx3 259 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1 260 fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)+dfw1 261 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2 262 fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)+dfw2 263 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3 264 fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)+dfw3 265 if(ip2(8)) 266 + ep2(1)=ep2(1)+half*ang(iwh,4,2)*(rww-ang(iwh,3,2))**2 267 if(ip3(8)) 268 + ep3(1)=ep3(1)+half*ang(iwh,4,3)*(rww-ang(iwh,3,3))**2 269 if(ith(8)) then 270 deriv(8,1)=deriv(8,1)+ 271 + dbond*(half*dbond*ang(iwh,4,4)-for*ang(iwh,3,4)) 272 endif 273 150 continue 274 rub(iwh,2)=rub(iwh,2)+eub 275 eww(13,1)=eww(13,1)+eub 276 if(ip2(8)) ep2(1)=ep2(1)-eub 277 if(ip3(8)) ep3(1)=ep3(1)-eub 278 1140 continue 279 endif 280 do 70 iwd=1,nwd 281 iwa=idih(iwd,1) 282 jwa=idih(iwd,2) 283 kwa=idih(iwd,3) 284 lwa=idih(iwd,4) 285 angle=dih(iwd,2,iset) 286 for=dih(iwd,3,iset) 287 rmul=dih(iwd,1,iset) 288 rdih(iwd,2)=zero 289 do 80 iwm=iwfr,iwto 290 xwij1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 291 xwij2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 292 xwij3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 293 xwkj1=xw(iwm,1,kwa)-xw(iwm,1,jwa) 294 xwkj2=xw(iwm,2,kwa)-xw(iwm,2,jwa) 295 xwkj3=xw(iwm,3,kwa)-xw(iwm,3,jwa) 296 xwkl1=xw(iwm,1,kwa)-xw(iwm,1,lwa) 297 xwkl2=xw(iwm,2,kwa)-xw(iwm,2,lwa) 298 xwkl3=xw(iwm,3,kwa)-xw(iwm,3,lwa) 299 xwik1=xwij1-xwkj1 300 xwik2=xwij2-xwkj2 301 xwik3=xwij3-xwkj3 302 xwjl1=xwkl1-xwkj1 303 xwjl2=xwkl2-xwkj2 304 xwjl3=xwkl3-xwkj3 305 xm1=xwij2*xwkj3-xwkj2*xwij3 306 xm2=xwij3*xwkj1-xwkj3*xwij1 307 xm3=xwij1*xwkj2-xwkj1*xwij2 308 xn1=xwkj2*xwkl3-xwkl2*xwkj3 309 xn2=xwkj3*xwkl1-xwkl3*xwkj1 310 xn3=xwkj1*xwkl2-xwkl1*xwkj2 311 rm2i=one/(xm1**2+xm2**2+xm3**2) 312 rn2i=one/(xn1**2+xn2**2+xn3**2) 313 rmni=sqrt(rm2i*rn2i) 314 cphi=(xm1*xn1+xm2*xn2+xm3*xn3)*rmni 315 if(cphi.lt.-one) cphi=-one 316 if(cphi.gt. one) cphi= one 317 phi=acos(cphi) 318 s=xwkj1*(xm2*xn3-xm3*xn2) +xwkj2*(xm3*xn1-xm1*xn3) 319 + +xwkj3*(xm1*xn2-xm2*xn1) 320 if(s.lt.zero) phi=-phi 321 sphi=sin(phi) 322 rpa=rmul*phi-angle 323 if(iand(iwdt(iwm),mdynam).eq.ldynam) 324 + rdih(iwd,2)=rdih(iwd,2)+for*(one+cos(rpa)) 325 dfor=(-for)*rmul*sin(rpa) 326 if(ip2(8)) ep2(1)=ep2(1)+ 327 + dih(iwd,3,2)*(one+cos(dih(iwd,1,2)*phi-dih(iwd,2,2))) 328 if(ip3(8)) ep3(1)=ep3(1)+ 329 + dih(iwd,3,3)*(one+cos(dih(iwd,1,3)*phi-dih(iwd,2,3))) 330 if(abs(sphi).lt.small) sphi=sign(small,sphi) 331 sphii=one/sphi 332 xd1=(-dfor)*sphii*(rmni*xn1-cphi*rm2i*xm1) 333 xe1=(-dfor)*sphii*(rmni*xm1-cphi*rn2i*xn1) 334 xd2=(-dfor)*sphii*(rmni*xn2-cphi*rm2i*xm2) 335 xe2=(-dfor)*sphii*(rmni*xm2-cphi*rn2i*xn2) 336 xd3=(-dfor)*sphii*(rmni*xn3-cphi*rm2i*xm3) 337 xe3=(-dfor)*sphii*(rmni*xm3-cphi*rn2i*xn3) 338 dfwi1=xwkj2*xd3-xwkj3*xd2 339 dfwi2=xwkj3*xd1-xwkj1*xd3 340 dfwi3=xwkj1*xd2-xwkj2*xd1 341 dfwj1=xwik2*xd3-xwik3*xd2-xwkl2*xe3+xwkl3*xe2 342 dfwj2=xwik3*xd1-xwik1*xd3-xwkl3*xe1+xwkl1*xe3 343 dfwj3=xwik1*xd2-xwik2*xd1-xwkl1*xe2+xwkl2*xe1 344 dfwk1=xwjl2*xe3-xwjl3*xe2-xwij2*xd3+xwij3*xd2 345 dfwk2=xwjl3*xe1-xwjl1*xe3-xwij3*xd1+xwij1*xd3 346 dfwk3=xwjl1*xe2-xwjl2*xe1-xwij1*xd2+xwij2*xd1 347 dfwl1=xwkj2*xe3-xwkj3*xe2 348 dfwl2=xwkj3*xe1-xwkj1*xe3 349 dfwl3=xwkj1*xe2-xwkj2*xe1 350 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfwi1 351 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfwi2 352 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfwi3 353 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfwj1 354 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfwj2 355 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfwj3 356 fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)-dfwk1 357 fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)-dfwk2 358 fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)-dfwk3 359 fw(iwm,1,lwa,1)=fw(iwm,1,lwa,1)-dfwl1 360 fw(iwm,2,lwa,1)=fw(iwm,2,lwa,1)-dfwl2 361 fw(iwm,3,lwa,1)=fw(iwm,3,lwa,1)-dfwl3 362 if(ith(9)) then 363 deriv(9,1)=deriv(9,1)+(one+cos(rpa))*dih(iwd,3,4) 364 + -for*sin(rpa)*(phi*dih(iwd,1,4)-dih(iwd,2,4)) 365 endif 366 80 continue 367 eww(3,1)=eww(3,1)+rdih(iwd,2) 368 if(ip2(8)) ep2(1)=ep2(1)-rdih(iwd,2) 369 if(ip3(8)) ep3(1)=ep3(1)-rdih(iwd,2) 370 70 continue 371 do 90 iwo=1,nwo 372 iwa=iimp(iwo,1) 373 jwa=iimp(iwo,2) 374 kwa=iimp(iwo,3) 375 lwa=iimp(iwo,4) 376 angle=dimp(iwo,2,iset) 377 for=dimp(iwo,3,iset) 378 rimp(iwo,2)=zero 379 do 100 iwm=iwfr,iwto 380 xwij1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 381 xwij2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 382 xwij3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 383 xwkj1=xw(iwm,1,kwa)-xw(iwm,1,jwa) 384 xwkj2=xw(iwm,2,kwa)-xw(iwm,2,jwa) 385 xwkj3=xw(iwm,3,kwa)-xw(iwm,3,jwa) 386 xwkl1=xw(iwm,1,kwa)-xw(iwm,1,lwa) 387 xwkl2=xw(iwm,2,kwa)-xw(iwm,2,lwa) 388 xwkl3=xw(iwm,3,kwa)-xw(iwm,3,lwa) 389 xwik1=xwij1-xwkj1 390 xwik2=xwij2-xwkj2 391 xwik3=xwij3-xwkj3 392 xwjl1=xwkl1-xwkj1 393 xwjl2=xwkl2-xwkj2 394 xwjl3=xwkl3-xwkj3 395 xm1=xwij2*xwkj3-xwkj2*xwij3 396 xm2=xwij3*xwkj1-xwkj3*xwij1 397 xm3=xwij1*xwkj2-xwkj1*xwij2 398 xn1=xwkj2*xwkl3-xwkl2*xwkj3 399 xn2=xwkj3*xwkl1-xwkl3*xwkj1 400 xn3=xwkj1*xwkl2-xwkl1*xwkj2 401 rm2i=one/(xm1**2+xm2**2+xm3**2) 402 rn2i=one/(xn1**2+xn2**2+xn3**2) 403 rmni=sqrt(rm2i*rn2i) 404 cphi=(xm1*xn1+xm2*xn2+xm3*xn3) 405 if(cphi.lt.-one) cphi=-one 406 if(cphi.gt. one) cphi= one 407 phi=acos(cphi) 408 s=xwkj1*(xm2*xn3-xm3*xn2) +xwkj2*(xm3*xn1-xm1*xn3) 409 + +xwkj3*(xm1*xn2-xm2*xn1) 410 if(s.lt.zero) phi=-phi 411 sphi=sin(phi) 412 dangle=(phi-angle)-nint((phi-angle)/twopi)*twopi 413 dfor=for*dangle 414 if(iand(iwdt(iwm),mdynam).eq.ldynam) rimp(iwo,2)=half*dfor*dangle 415 if(ip2(9)) then 416 danglep=(phi-dimp(iwo,2,2))-nint((phi-dimp(iwo,2,2))/twopi)*twopi 417 ep2(1)=ep2(1)+half*dimp(iwo,3,2)*danglep**2 418 endif 419 if(ip3(9)) then 420 danglep=(phi-dimp(iwo,2,3))-nint((phi-dimp(iwo,2,3))/twopi)*twopi 421 ep3(1)=ep3(1)+half*dimp(iwo,3,3)*danglep**2 422 endif 423 if(abs(sphi).lt.small) sphi=sign(small,sphi) 424 sphii=one/sphi 425 xd1=(-dfor)*sphii*(rmni*xn1-cphi*rm2i*xm1) 426 xe1=(-dfor)*sphii*(rmni*xm1-cphi*rn2i*xn1) 427 xd2=(-dfor)*sphii*(rmni*xn2-cphi*rm2i*xm2) 428 xe2=(-dfor)*sphii*(rmni*xm2-cphi*rn2i*xn2) 429 xd3=(-dfor)*sphii*(rmni*xn3-cphi*rm2i*xm3) 430 xe3=(-dfor)*sphii*(rmni*xm3-cphi*rn2i*xn3) 431 dfwi1=xwkj2*xd3-xwkj3*xd2 432 dfwi2=xwkj3*xd1-xwkj1*xd3 433 dfwi3=xwkj1*xd2-xwkj2*xd1 434 dfwj1=xwik2*xd3-xwik3*xd2-xwkl2*xe3+xwkl3*xe2 435 dfwj2=xwik3*xd1-xwik1*xd3-xwkl3*xe1+xwkl1*xe3 436 dfwj3=xwik1*xd2-xwik2*xd1-xwkl1*xe2+xwkl2*xe1 437 dfwk1=xwjl2*xe3-xwjl3*xe2-xwij2*xd3+xwij3*xd2 438 dfwk2=xwjl3*xe1-xwjl1*xe3-xwij3*xd1+xwij1*xd3 439 dfwk3=xwjl1*xe2-xwjl2*xe1-xwij1*xd2+xwij2*xd1 440 dfwl1=xwkj2*xe3-xwkj3*xe2 441 dfwl2=xwkj3*xe1-xwkj1*xe3 442 dfwl3=xwkj1*xe2-xwkj2*xe1 443 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfwi1 444 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfwi2 445 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfwi3 446 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfwj1 447 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfwj2 448 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfwj3 449 fw(iwm,1,kwa,1)=fw(iwm,1,kwa,1)-dfwk1 450 fw(iwm,2,kwa,1)=fw(iwm,2,kwa,1)-dfwk2 451 fw(iwm,3,kwa,1)=fw(iwm,3,kwa,1)-dfwk3 452 fw(iwm,1,lwa,1)=fw(iwm,1,lwa,1)-dfwl1 453 fw(iwm,2,lwa,1)=fw(iwm,2,lwa,1)-dfwl2 454 fw(iwm,3,lwa,1)=fw(iwm,3,lwa,1)-dfwl3 455 if(ith(10)) then 456 deriv(10,1)=deriv(10,1)+ 457 + dangle*(half*dangle*dimp(iwo,3,4)-for*dimp(iwo,2,4)) 458 endif 459 100 continue 460 eww(4,1)=eww(4,1)+rimp(iwo,2) 461 if(ip2(9)) ep2(1)=ep2(1)-rimp(iwo,2) 462 if(ip3(9)) ep3(1)=ep3(1)-rimp(iwo,2) 463 90 continue 464 c6p1=zero 465 c12p1=zero 466 c6p2=zero 467 c12p2=zero 468 qip1=zero 469 qjp1=zero 470 qip2=zero 471 qjp2=zero 472 do 110 iwt=1,nwt 473 iwa=idwt(iwt,1) 474 jwa=idwt(iwt,2) 475 c6=vdw(iwatm(iwa),iwatm(jwa),2,iset) 476 c12=vdw(iwatm(iwa),iwatm(jwa),4,iset) 477 if(ip2(2)) then 478 c6p1=vdw(iwatm(iwa),iwatm(jwa),2,2) 479 c12p1=vdw(iwatm(iwa),iwatm(jwa),4,2) 480 endif 481 if(ip3(2)) then 482 c6p2=vdw(iwatm(iwa),iwatm(jwa),2,3) 483 c12p2=vdw(iwatm(iwa),iwatm(jwa),4,3) 484 endif 485 if(ith(2).or.ith(4)) then 486 c6t=vdw(iwatm(iwa),iwatm(jwa),2,4) 487 c12t=vdw(iwatm(iwa),iwatm(jwa),4,4) 488 qit=chg(iwq(iwa),1,4)*q14fac 489 qjt=chg(iwq(jwa),1,4) 490 endif 491 cf6=six*c6 492 cf12=twelve*c12 493 qi=chg(iwq(iwa),1,iset)*q14fac 494 qj=chg(iwq(jwa),1,iset) 495 q=qi*qj 496 if(ip2(4)) then 497 qip1=chg(iwq(iwa),1,2)*q14fac 498 qjp1=chg(iwq(jwa),1,2) 499 qp1=qip1*qjp1 500 endif 501 if(ip3(4)) then 502 qip2=chg(iwq(iwa),1,3)*q14fac 503 qjp2=chg(iwq(jwa),1,3) 504 qp2=qip2*qjp2 505 endif 506 ep2l=zero 507 ep3l=zero 508 ep2q=zero 509 ep3q=zero 510 do 120 iwm=iwfr,iwto 511 rxx=xw(iwm,1,iwa)-xw(iwm,1,jwa) 512 rxy=xw(iwm,2,iwa)-xw(iwm,2,jwa) 513 rxz=xw(iwm,3,iwa)-xw(iwm,3,jwa) 514 r2=rxx*rxx+rxy*rxy+rxz*rxz 515 r2i=one/r2 516 r1i=sqrt(r2i) 517 r6i=r2i*r2i*r2i 518 eterml=(c12*r6i-c6)*r6i 519 etermq=q*r1i 520 if(iand(iwdt(iwm),mdynam).eq.ldynam) eww(5,1)=eww(5,1)+eterml 521 if(iand(iwdt(iwm),mdynam).eq.ldynam) eww(6,1)=eww(6,1)+etermq 522 if(ip2(2)) ep2l=ep2l-eterml+(c12p1*r6i-c6p1)*r6i 523 if(ip3(2)) ep3l=ep3l-eterml+(c12p2*r6i-c6p2)*r6i 524 if(ip2(4)) ep2q=ep2q-etermq+qp1*r1i 525 if(ip3(4)) ep3q=ep3q-etermq+qp2*r1i 526 dfw=((cf12*r6i-cf6)*r6i+q*r1i)*r2i 527 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)+dfw*rxx 528 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)+dfw*rxy 529 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)+dfw*rxz 530 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw*rxx 531 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw*rxy 532 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw*rxz 533 if(ith(2)) then 534 deriv(2,1)=deriv(2,1)+(c12t*r6i-c6t)*r6i 535 endif 536 if(ith(4)) then 537 deriv(4,1)=deriv(4,1)+(qi*qjt+qj*qit)*r1i 538 endif 539 120 continue 540 ep2(1)=ep2(1)+ep2l+ep2q 541 ep3(1)=ep3(1)+ep3l+ep3q 542 if(ipme.ne.0) then 543 qij=(one-q14fac)*chg(iwq(iwa),1,iset)*chg(iwq(jwa),1,iset) 544 do 111 iwm=iwfr,iwto 545 rwx1=xw(iwm,1,iwa)-xw(iwm,1,jwa) 546 rwx2=xw(iwm,2,iwa)-xw(iwm,2,jwa) 547 rwx3=xw(iwm,3,iwa)-xw(iwm,3,jwa) 548 rww=sqrt(rwx1**2+rwx2**2+rwx3**2) 549 rwi=one/rww 550 ferfc=one-erfc(ealpha*rww) 551 fderfc=-(ealpha*derfc(ealpha*rww)) 552 epmecw=epmecw-ferfc*qij*rwi 553 eww(6,1)=eww(6,1)-ferfc*qij*rwi 554 dfor=-(qij*rwi*rwi*(ferfc*rwi-fderfc)) 555 dfw1=dfor*rwx1 556 dfw2=dfor*rwx2 557 dfw3=dfor*rwx3 558 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)-dfw1 559 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)+dfw1 560 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)-dfw2 561 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)+dfw2 562 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)-dfw3 563 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)+dfw3 564 vpmeb(1)=vpmeb(1)+dfw1*rwx1 565 vpmeb(2)=vpmeb(2)+dfw2*rwx1 566 vpmeb(3)=vpmeb(3)+dfw3*rwx1 567 vpmeb(4)=vpmeb(4)+dfw2*rwx2 568 vpmeb(5)=vpmeb(5)+dfw3*rwx2 569 vpmeb(6)=vpmeb(6)+dfw3*rwx3 570 111 continue 571 endif 572 110 continue 573 do 130 iwn=1,nwn 574 iwa=idwn(iwn,1) 575 jwa=idwn(iwn,2) 576 c6=vdw(iwatm(iwa),iwatm(jwa),1,iset) 577 c12=vdw(iwatm(iwa),iwatm(jwa),3,iset) 578 if(ip2(2)) then 579 c6p1=vdw(iwatm(iwa),iwatm(jwa),1,2) 580 c12p1=vdw(iwatm(iwa),iwatm(jwa),3,2) 581 endif 582 if(ip3(2)) then 583 c6p2=vdw(iwatm(iwa),iwatm(jwa),1,3) 584 c12p2=vdw(iwatm(iwa),iwatm(jwa),3,3) 585 endif 586 if(ith(2).or.ith(4)) then 587 c6t=vdw(iwatm(iwa),iwatm(jwa),1,4) 588 c12t=vdw(iwatm(iwa),iwatm(jwa),3,4) 589 qit=chg(iwq(iwa),1,4) 590 qjt=chg(iwq(jwa),1,4) 591 endif 592 cf6=six*c6 593 cf12=twelve*c12 594 qi=chg(iwq(iwa),1,iset) 595 qj=chg(iwq(jwa),1,iset) 596 q=qi*qj 597 if(ip2(4)) then 598 qip1=chg(iwq(iwa),1,2) 599 qjp1=chg(iwq(jwa),1,2) 600 qp1=qip1*qjp1 601 endif 602 if(ip3(4)) then 603 qip2=chg(iwq(iwa),1,3) 604 qjp2=chg(iwq(jwa),1,3) 605 qp2=qip2*qjp2 606 endif 607 ep2l=zero 608 ep3l=zero 609 ep2q=zero 610 ep3q=zero 611 do 140 iwm=iwfr,iwto 612 rxx=xw(iwm,1,iwa)-xw(iwm,1,jwa) 613 rxy=xw(iwm,2,iwa)-xw(iwm,2,jwa) 614 rxz=xw(iwm,3,iwa)-xw(iwm,3,jwa) 615 r2=rxx*rxx+rxy*rxy+rxz*rxz 616 r2i=one/r2 617 r1i=sqrt(r2i) 618 r6i=r2i*r2i*r2i 619 ferfc=one 620 fderfc=zero 621 if(ipme.ne.0) then 622 ferfc=erfc(ealpha/r1i) 623 fderfc=ealpha+derfc(ealpha/r1i) 624 endif 625 eterml=(c12*r6i-c6)*r6i 626 etermq=ferfc*q*r1i 627 if(iand(iwdt(iwm),mdynam).eq.ldynam) then 628 eww(5,1)=eww(5,1)+eterml 629 eww(6,1)=eww(6,1)+etermq 630 endif 631 if(ip2(2)) ep2l=ep2l-eterml+(c12p1*r6i-c6p1)*r6i 632 if(ip3(2)) ep3l=ep3l-eterml+(c12p2*r6i-c6p2)*r6i 633 if(ip2(4)) ep2q=ep2q-etermq+qp1*r1i 634 if(ip3(4)) ep3q=ep3q-etermq+qp2*r1i 635 dfw=((cf12*r6i-cf6)*r6i+q*(ferfc*r1i-fderfc))*r2i 636 fw(iwm,1,iwa,1)=fw(iwm,1,iwa,1)+dfw*rxx 637 fw(iwm,2,iwa,1)=fw(iwm,2,iwa,1)+dfw*rxy 638 fw(iwm,3,iwa,1)=fw(iwm,3,iwa,1)+dfw*rxz 639 fw(iwm,1,jwa,1)=fw(iwm,1,jwa,1)-dfw*rxx 640 fw(iwm,2,jwa,1)=fw(iwm,2,jwa,1)-dfw*rxy 641 fw(iwm,3,jwa,1)=fw(iwm,3,jwa,1)-dfw*rxz 642 if(ith(2)) deriv(2,1)=deriv(2,1)+(c12t*r6i-c6t)*r6i 643 if(ith(4)) deriv(4,1)=deriv(4,1)+(qi*qjt+qj*qit)*r1i 644 140 continue 645 ep2(1)=ep2(1)+ep2l+ep2q 646 ep3(1)=ep3(1)+ep3l+ep3q 647 130 continue 648c 649 return 650 end 651