1c Extern "C" declaration has the form: 2c 3c void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *, 4c int *, int *, int *, int *, double *, double *, 5c double *, double *, double *, double *, double *, double *, 6c double *, double *, double *, double *, double *, double *, 7c double *, double *, double *, double *, double *, double *, int *); 8c 9c Call from pair_meam.cpp has the form: 10c 11c meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom, 12c &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], 13c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], 14c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], 15c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, 16c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], 17c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); 18c 19 20 subroutine meam_force(i, nmax, 21 $ eflag_either, eflag_global, eflag_atom, vflag_atom, 22 $ eng_vdwl, eatom, ntype, type, fmap, x, 23 $ numneigh, firstneigh, numneigh_full, firstneigh_full, 24 $ scrfcn, dscrfcn, fcpair, 25 $ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp, 26 $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f, 27 $ vatom, errorflag) 28 29 use meam_data 30 implicit none 31 32 integer eflag_either, eflag_global, eflag_atom, vflag_atom 33 integer nmax, ntype, type, fmap 34 real*8 eng_vdwl, eatom, x 35 integer numneigh, firstneigh, numneigh_full, firstneigh_full 36 real*8 scrfcn, dscrfcn, fcpair 37 real*8 dGamma1, dGamma2, dGamma3 38 real*8 rho0, rho1, rho2, rho3, fp 39 real*8 Arho1, Arho2, Arho2b 40 real*8 Arho3, Arho3b 41 real*8 t_ave, tsq_ave, f, vatom 42 integer errorflag 43 44 dimension eatom(nmax) 45 dimension type(nmax), fmap(ntype) 46 dimension x(3,nmax) 47 dimension firstneigh(numneigh), firstneigh_full(numneigh_full) 48 dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh) 49 dimension dGamma1(nmax), dGamma2(nmax), dGamma3(nmax) 50 dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax) 51 dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) 52 dimension Arho3(10,nmax), Arho3b(3,nmax) 53 dimension t_ave(3,nmax), tsq_ave(3,nmax), f(3,nmax), vatom(6,nmax) 54 55 integer i,j,jn,k,kn,kk,m,n,p,q 56 integer nv2,nv3,elti,eltj,eltk,ind 57 real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3 58 real*8 delik(3),deljk(3),v(6),fi(3),fj(3) 59 real*8 Eu,astar,astarp,third,sixth 60 real*8 pp,phiforce,dUdrij,dUdsij,dUdrijm(3),force,forcem 61 real*8 B,r,recip,phi,phip,rhop,a 62 real*8 sij,fcij,dfcij,ds(3) 63 real*8 a0,a1,a1i,a1j,a2,a2i,a2j 64 real*8 a3i,a3j,a3i1,a3i2,a3j1,a3j2 65 real*8 G,dG,Gbar,dGbar,gam,shpi(3),shpj(3),Z,denom 66 real*8 ai,aj,ro0i,ro0j,invrei,invrej 67 real*8 b0,rhoa0j,drhoa0j,rhoa0i,drhoa0i 68 real*8 b1,rhoa1j,drhoa1j,rhoa1i,drhoa1i 69 real*8 b2,rhoa2j,drhoa2j,rhoa2i,drhoa2i 70 real*8 a3,a3a,b3,rhoa3j,drhoa3j,rhoa3i,drhoa3i 71 real*8 drho0dr1,drho0dr2,drho0ds1,drho0ds2 72 real*8 drho1dr1,drho1dr2,drho1ds1,drho1ds2 73 real*8 drho1drm1(3),drho1drm2(3) 74 real*8 drho2dr1,drho2dr2,drho2ds1,drho2ds2 75 real*8 drho2drm1(3),drho2drm2(3) 76 real*8 drho3dr1,drho3dr2,drho3ds1,drho3ds2 77 real*8 drho3drm1(3),drho3drm2(3) 78 real*8 dt1dr1,dt1dr2,dt1ds1,dt1ds2 79 real*8 dt2dr1,dt2dr2,dt2ds1,dt2ds2 80 real*8 dt3dr1,dt3dr2,dt3ds1,dt3ds2 81 real*8 drhodr1,drhodr2,drhods1,drhods2,drhodrm1(3),drhodrm2(3) 82 real*8 arg,arg1,arg2 83 real*8 arg1i1,arg1j1,arg1i2,arg1j2,arg2i2,arg2j2 84 real*8 arg1i3,arg1j3,arg2i3,arg2j3,arg3i3,arg3j3 85 real*8 dsij1,dsij2,force1,force2 86 real*8 t1i,t2i,t3i,t1j,t2j,t3j 87 88 errorflag = 0 89 third = 1.0/3.0 90 sixth = 1.0/6.0 91 92c Compute forces atom i 93 94 elti = fmap(type(i)) 95 96 if (elti.gt.0) then 97 xitmp = x(1,i) 98 yitmp = x(2,i) 99 zitmp = x(3,i) 100 101c Treat each pair 102 do jn = 1,numneigh 103 104 j = firstneigh(jn) 105 eltj = fmap(type(j)) 106 107 if (scrfcn(jn).ne.0.d0.and.eltj.gt.0) then 108 109 sij = scrfcn(jn)*fcpair(jn) 110 delij(1) = x(1,j) - xitmp 111 delij(2) = x(2,j) - yitmp 112 delij(3) = x(3,j) - zitmp 113 rij2 = delij(1)*delij(1) + delij(2)*delij(2) 114 $ + delij(3)*delij(3) 115 if (rij2.lt.cutforcesq) then 116 rij = sqrt(rij2) 117 r = rij 118 119c Compute phi and phip 120 ind = eltind(elti,eltj) 121 pp = rij*rdrar + 1.0D0 122 kk = pp 123 kk = min(kk,nrar-1) 124 pp = pp - kk 125 pp = min(pp,1.0D0) 126 phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp 127 $ + phirar1(kk,ind))*pp + phirar(kk,ind) 128 phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp 129 $ + phirar4(kk,ind) 130 recip = 1.0d0/r 131 132 if (eflag_either.ne.0) then 133 if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij 134 if (eflag_atom.ne.0) then 135 eatom(i) = eatom(i) + 0.5*phi*sij 136 eatom(j) = eatom(j) + 0.5*phi*sij 137 endif 138 endif 139 140c write(1,*) "force_meamf: phi: ",phi 141c write(1,*) "force_meamf: phip: ",phip 142 143c Compute pair densities and derivatives 144 invrei = 1.d0/re_meam(elti,elti) 145 ai = rij*invrei - 1.d0 146 ro0i = rho0_meam(elti) 147 rhoa0i = ro0i*exp(-beta0_meam(elti)*ai) 148 drhoa0i = -beta0_meam(elti)*invrei*rhoa0i 149 rhoa1i = ro0i*exp(-beta1_meam(elti)*ai) 150 drhoa1i = -beta1_meam(elti)*invrei*rhoa1i 151 rhoa2i = ro0i*exp(-beta2_meam(elti)*ai) 152 drhoa2i = -beta2_meam(elti)*invrei*rhoa2i 153 rhoa3i = ro0i*exp(-beta3_meam(elti)*ai) 154 drhoa3i = -beta3_meam(elti)*invrei*rhoa3i 155 156 if (elti.ne.eltj) then 157 invrej = 1.d0/re_meam(eltj,eltj) 158 aj = rij*invrej - 1.d0 159 ro0j = rho0_meam(eltj) 160 rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj) 161 drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j 162 rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj) 163 drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j 164 rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj) 165 drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j 166 rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj) 167 drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j 168 else 169 rhoa0j = rhoa0i 170 drhoa0j = drhoa0i 171 rhoa1j = rhoa1i 172 drhoa1j = drhoa1i 173 rhoa2j = rhoa2i 174 drhoa2j = drhoa2i 175 rhoa3j = rhoa3i 176 drhoa3j = drhoa3i 177 endif 178 179 if (ialloy.eq.1) then 180 rhoa1j = rhoa1j * t1_meam(eltj) 181 rhoa2j = rhoa2j * t2_meam(eltj) 182 rhoa3j = rhoa3j * t3_meam(eltj) 183 rhoa1i = rhoa1i * t1_meam(elti) 184 rhoa2i = rhoa2i * t2_meam(elti) 185 rhoa3i = rhoa3i * t3_meam(elti) 186 drhoa1j = drhoa1j * t1_meam(eltj) 187 drhoa2j = drhoa2j * t2_meam(eltj) 188 drhoa3j = drhoa3j * t3_meam(eltj) 189 drhoa1i = drhoa1i * t1_meam(elti) 190 drhoa2i = drhoa2i * t2_meam(elti) 191 drhoa3i = drhoa3i * t3_meam(elti) 192 endif 193 194 nv2 = 1 195 nv3 = 1 196 arg1i1 = 0.d0 197 arg1j1 = 0.d0 198 arg1i2 = 0.d0 199 arg1j2 = 0.d0 200 arg1i3 = 0.d0 201 arg1j3 = 0.d0 202 arg3i3 = 0.d0 203 arg3j3 = 0.d0 204 do n = 1,3 205 do p = n,3 206 do q = p,3 207 arg = delij(n)*delij(p)*delij(q)*v3D(nv3) 208 arg1i3 = arg1i3 + Arho3(nv3,i)*arg 209 arg1j3 = arg1j3 - Arho3(nv3,j)*arg 210 nv3 = nv3+1 211 enddo 212 arg = delij(n)*delij(p)*v2D(nv2) 213 arg1i2 = arg1i2 + Arho2(nv2,i)*arg 214 arg1j2 = arg1j2 + Arho2(nv2,j)*arg 215 nv2 = nv2+1 216 enddo 217 arg1i1 = arg1i1 + Arho1(n,i)*delij(n) 218 arg1j1 = arg1j1 - Arho1(n,j)*delij(n) 219 arg3i3 = arg3i3 + Arho3b(n,i)*delij(n) 220 arg3j3 = arg3j3 - Arho3b(n,j)*delij(n) 221 enddo 222 223c rho0 terms 224 drho0dr1 = drhoa0j * sij 225 drho0dr2 = drhoa0i * sij 226 227c rho1 terms 228 a1 = 2*sij/rij 229 drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1 230 drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1 231 a1 = 2.d0*sij/rij 232 do m = 1,3 233 drho1drm1(m) = a1*rhoa1j*Arho1(m,i) 234 drho1drm2(m) = -a1*rhoa1i*Arho1(m,j) 235 enddo 236 237c rho2 terms 238 a2 = 2*sij/rij2 239 drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2 240 $ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij 241 drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2 242 $ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij 243 a2 = 4*sij/rij2 244 do m = 1,3 245 drho2drm1(m) = 0.d0 246 drho2drm2(m) = 0.d0 247 do n = 1,3 248 drho2drm1(m) = drho2drm1(m) 249 $ + Arho2(vind2D(m,n),i)*delij(n) 250 drho2drm2(m) = drho2drm2(m) 251 $ - Arho2(vind2D(m,n),j)*delij(n) 252 enddo 253 drho2drm1(m) = a2*rhoa2j*drho2drm1(m) 254 drho2drm2(m) = -a2*rhoa2i*drho2drm2(m) 255 enddo 256 257c rho3 terms 258 rij3 = rij*rij2 259 a3 = 2*sij/rij3 260 a3a = 6.d0/5.d0*sij/rij 261 drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3 262 $ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3 263 drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3 264 $ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3 265 a3 = 6*sij/rij3 266 a3a = 6*sij/(5*rij) 267 do m = 1,3 268 drho3drm1(m) = 0.d0 269 drho3drm2(m) = 0.d0 270 nv2 = 1 271 do n = 1,3 272 do p = n,3 273 arg = delij(n)*delij(p)*v2D(nv2) 274 drho3drm1(m) = drho3drm1(m) 275 $ + Arho3(vind3D(m,n,p),i)*arg 276 drho3drm2(m) = drho3drm2(m) 277 $ + Arho3(vind3D(m,n,p),j)*arg 278 nv2 = nv2 + 1 279 enddo 280 enddo 281 drho3drm1(m) = (a3*drho3drm1(m) - a3a*Arho3b(m,i)) 282 $ *rhoa3j 283 drho3drm2(m) = (-a3*drho3drm2(m) + a3a*Arho3b(m,j)) 284 $ *rhoa3i 285 enddo 286 287c Compute derivatives of weighting functions t wrt rij 288 t1i = t_ave(1,i) 289 t2i = t_ave(2,i) 290 t3i = t_ave(3,i) 291 t1j = t_ave(1,j) 292 t2j = t_ave(2,j) 293 t3j = t_ave(3,j) 294 295 if (ialloy.eq.1) then 296 297 a1i = 0.d0 298 a1j = 0.d0 299 a2i = 0.d0 300 a2j = 0.d0 301 a3i = 0.d0 302 a3j = 0.d0 303 if ( tsq_ave(1,i) .ne. 0.d0 ) then 304 a1i = drhoa0j*sij/tsq_ave(1,i) 305 endif 306 if ( tsq_ave(1,j) .ne. 0.d0 ) then 307 a1j = drhoa0i*sij/tsq_ave(1,j) 308 endif 309 if ( tsq_ave(2,i) .ne. 0.d0 ) then 310 a2i = drhoa0j*sij/tsq_ave(2,i) 311 endif 312 if ( tsq_ave(2,j) .ne. 0.d0 ) then 313 a2j = drhoa0i*sij/tsq_ave(2,j) 314 endif 315 if ( tsq_ave(3,i) .ne. 0.d0 ) then 316 a3i = drhoa0j*sij/tsq_ave(3,i) 317 endif 318 if ( tsq_ave(3,j) .ne. 0.d0 ) then 319 a3j = drhoa0i*sij/tsq_ave(3,j) 320 endif 321 322 dt1dr1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2) 323 dt1dr2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2) 324 dt2dr1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2) 325 dt2dr2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2) 326 dt3dr1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2) 327 dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) 328 329 else if (ialloy.eq.2) then 330 331 dt1dr1 = 0.d0 332 dt1dr2 = 0.d0 333 dt2dr1 = 0.d0 334 dt2dr2 = 0.d0 335 dt3dr1 = 0.d0 336 dt3dr2 = 0.d0 337 338 else 339 340 ai = 0.d0 341 if( rho0(i) .ne. 0.d0 ) then 342 ai = drhoa0j*sij/rho0(i) 343 end if 344 aj = 0.d0 345 if( rho0(j) .ne. 0.d0 ) then 346 aj = drhoa0i*sij/rho0(j) 347 end if 348 349 dt1dr1 = ai*(t1_meam(eltj)-t1i) 350 dt1dr2 = aj*(t1_meam(elti)-t1j) 351 dt2dr1 = ai*(t2_meam(eltj)-t2i) 352 dt2dr2 = aj*(t2_meam(elti)-t2j) 353 dt3dr1 = ai*(t3_meam(eltj)-t3i) 354 dt3dr2 = aj*(t3_meam(elti)-t3j) 355 356 endif 357 358c Compute derivatives of total density wrt rij, sij and rij(3) 359 call get_shpfcn(shpi,lattce_meam(elti,elti)) 360 call get_shpfcn(shpj,lattce_meam(eltj,eltj)) 361 drhodr1 = dGamma1(i)*drho0dr1 362 $ + dGamma2(i)* 363 $ (dt1dr1*rho1(i)+t1i*drho1dr1 364 $ + dt2dr1*rho2(i)+t2i*drho2dr1 365 $ + dt3dr1*rho3(i)+t3i*drho3dr1) 366 $ - dGamma3(i)* 367 $ (shpi(1)*dt1dr1+shpi(2)*dt2dr1+shpi(3)*dt3dr1) 368 drhodr2 = dGamma1(j)*drho0dr2 369 $ + dGamma2(j)* 370 $ (dt1dr2*rho1(j)+t1j*drho1dr2 371 $ + dt2dr2*rho2(j)+t2j*drho2dr2 372 $ + dt3dr2*rho3(j)+t3j*drho3dr2) 373 $ - dGamma3(j)* 374 $ (shpj(1)*dt1dr2+shpj(2)*dt2dr2+shpj(3)*dt3dr2) 375 do m = 1,3 376 drhodrm1(m) = 0.d0 377 drhodrm2(m) = 0.d0 378 drhodrm1(m) = dGamma2(i)* 379 $ (t1i*drho1drm1(m) 380 $ + t2i*drho2drm1(m) 381 $ + t3i*drho3drm1(m)) 382 drhodrm2(m) = dGamma2(j)* 383 $ (t1j*drho1drm2(m) 384 $ + t2j*drho2drm2(m) 385 $ + t3j*drho3drm2(m)) 386 enddo 387 388c Compute derivatives wrt sij, but only if necessary 389 if (dscrfcn(jn).ne.0.d0) then 390 drho0ds1 = rhoa0j 391 drho0ds2 = rhoa0i 392 a1 = 2.d0/rij 393 drho1ds1 = a1*rhoa1j*arg1i1 394 drho1ds2 = a1*rhoa1i*arg1j1 395 a2 = 2.d0/rij2 396 drho2ds1 = a2*rhoa2j*arg1i2 397 $ - 2.d0/3.d0*Arho2b(i)*rhoa2j 398 drho2ds2 = a2*rhoa2i*arg1j2 399 $ - 2.d0/3.d0*Arho2b(j)*rhoa2i 400 a3 = 2.d0/rij3 401 a3a = 6.d0/(5.d0*rij) 402 drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3 403 drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3 404 405 if (ialloy.eq.1) then 406 407 a1i = 0.d0 408 a1j = 0.d0 409 a2i = 0.d0 410 a2j = 0.d0 411 a3i = 0.d0 412 a3j = 0.d0 413 if ( tsq_ave(1,i) .ne. 0.d0 ) then 414 a1i = rhoa0j/tsq_ave(1,i) 415 endif 416 if ( tsq_ave(1,j) .ne. 0.d0 ) then 417 a1j = rhoa0i/tsq_ave(1,j) 418 endif 419 if ( tsq_ave(2,i) .ne. 0.d0 ) then 420 a2i = rhoa0j/tsq_ave(2,i) 421 endif 422 if ( tsq_ave(2,j) .ne. 0.d0 ) then 423 a2j = rhoa0i/tsq_ave(2,j) 424 endif 425 if ( tsq_ave(3,i) .ne. 0.d0 ) then 426 a3i = rhoa0j/tsq_ave(3,i) 427 endif 428 if ( tsq_ave(3,j) .ne. 0.d0 ) then 429 a3j = rhoa0i/tsq_ave(3,j) 430 endif 431 432 dt1ds1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2) 433 dt1ds2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2) 434 dt2ds1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2) 435 dt2ds2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2) 436 dt3ds1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2) 437 dt3ds2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) 438 439 else if (ialloy.eq.2) then 440 441 dt1ds1 = 0.d0 442 dt1ds2 = 0.d0 443 dt2ds1 = 0.d0 444 dt2ds2 = 0.d0 445 dt3ds1 = 0.d0 446 dt3ds2 = 0.d0 447 448 else 449 450 ai = 0.d0 451 if( rho0(i) .ne. 0.d0 ) then 452 ai = rhoa0j/rho0(i) 453 end if 454 aj = 0.d0 455 if( rho0(j) .ne. 0.d0 ) then 456 aj = rhoa0i/rho0(j) 457 end if 458 459 dt1ds1 = ai*(t1_meam(eltj)-t1i) 460 dt1ds2 = aj*(t1_meam(elti)-t1j) 461 dt2ds1 = ai*(t2_meam(eltj)-t2i) 462 dt2ds2 = aj*(t2_meam(elti)-t2j) 463 dt3ds1 = ai*(t3_meam(eltj)-t3i) 464 dt3ds2 = aj*(t3_meam(elti)-t3j) 465 466 endif 467 468 drhods1 = dGamma1(i)*drho0ds1 469 $ + dGamma2(i)* 470 $ (dt1ds1*rho1(i)+t1i*drho1ds1 471 $ + dt2ds1*rho2(i)+t2i*drho2ds1 472 $ + dt3ds1*rho3(i)+t3i*drho3ds1) 473 $ - dGamma3(i)* 474 $ (shpi(1)*dt1ds1+shpi(2)*dt2ds1+shpi(3)*dt3ds1) 475 drhods2 = dGamma1(j)*drho0ds2 476 $ + dGamma2(j)* 477 $ (dt1ds2*rho1(j)+t1j*drho1ds2 478 $ + dt2ds2*rho2(j)+t2j*drho2ds2 479 $ + dt3ds2*rho3(j)+t3j*drho3ds2) 480 $ - dGamma3(j)* 481 $ (shpj(1)*dt1ds2+shpj(2)*dt2ds2+shpj(3)*dt3ds2) 482 endif 483 484c Compute derivatives of energy wrt rij, sij and rij(3) 485 dUdrij = phip*sij 486 $ + fp(i)*drhodr1 + fp(j)*drhodr2 487 dUdsij = 0.d0 488 if (dscrfcn(jn).ne.0.d0) then 489 dUdsij = phi 490 $ + fp(i)*drhods1 + fp(j)*drhods2 491 endif 492 do m = 1,3 493 dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m) 494 enddo 495 496c Add the part of the force due to dUdrij and dUdsij 497 498 force = dUdrij*recip + dUdsij*dscrfcn(jn) 499 do m = 1,3 500 forcem = delij(m)*force + dUdrijm(m) 501 f(m,i) = f(m,i) + forcem 502 f(m,j) = f(m,j) - forcem 503 enddo 504 505c Tabulate per-atom virial as symmetrized stress tensor 506 507 if (vflag_atom.ne.0) then 508 fi(1) = delij(1)*force + dUdrijm(1) 509 fi(2) = delij(2)*force + dUdrijm(2) 510 fi(3) = delij(3)*force + dUdrijm(3) 511 v(1) = -0.5 * (delij(1) * fi(1)) 512 v(2) = -0.5 * (delij(2) * fi(2)) 513 v(3) = -0.5 * (delij(3) * fi(3)) 514 v(4) = -0.25 * (delij(1)*fi(2) + delij(2)*fi(1)) 515 v(5) = -0.25 * (delij(1)*fi(3) + delij(3)*fi(1)) 516 v(6) = -0.25 * (delij(2)*fi(3) + delij(3)*fi(2)) 517 518 vatom(1,i) = vatom(1,i) + v(1) 519 vatom(2,i) = vatom(2,i) + v(2) 520 vatom(3,i) = vatom(3,i) + v(3) 521 vatom(4,i) = vatom(4,i) + v(4) 522 vatom(5,i) = vatom(5,i) + v(5) 523 vatom(6,i) = vatom(6,i) + v(6) 524 vatom(1,j) = vatom(1,j) + v(1) 525 vatom(2,j) = vatom(2,j) + v(2) 526 vatom(3,j) = vatom(3,j) + v(3) 527 vatom(4,j) = vatom(4,j) + v(4) 528 vatom(5,j) = vatom(5,j) + v(5) 529 vatom(6,j) = vatom(6,j) + v(6) 530 endif 531 532c Now compute forces on other atoms k due to change in sij 533 534 if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100 535 do kn = 1,numneigh_full 536 k = firstneigh_full(kn) 537 eltk = fmap(type(k)) 538 if (k.ne.j.and.eltk.gt.0) then 539 call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2, 540 $ ntype,type,fmap,x,scrfcn,fcpair) 541 if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then 542 force1 = dUdsij*dsij1 543 force2 = dUdsij*dsij2 544 do m = 1,3 545 delik(m) = x(m,k) - x(m,i) 546 deljk(m) = x(m,k) - x(m,j) 547 enddo 548 do m = 1,3 549 f(m,i) = f(m,i) + force1*delik(m) 550 f(m,j) = f(m,j) + force2*deljk(m) 551 f(m,k) = f(m,k) - force1*delik(m) 552 $ - force2*deljk(m) 553 enddo 554 555c Tabulate per-atom virial as symmetrized stress tensor 556 557 if (vflag_atom.ne.0) then 558 fi(1) = force1*delik(1) 559 fi(2) = force1*delik(2) 560 fi(3) = force1*delik(3) 561 fj(1) = force2*deljk(1) 562 fj(2) = force2*deljk(2) 563 fj(3) = force2*deljk(3) 564 v(1) = -third * (delik(1)*fi(1) + deljk(1)*fj(1)) 565 v(2) = -third * (delik(2)*fi(2) + deljk(2)*fj(2)) 566 v(3) = -third * (delik(3)*fi(3) + deljk(3)*fj(3)) 567 v(4) = -sixth * (delik(1)*fi(2) + deljk(1)*fj(2) + 568 $ delik(2)*fi(1) + deljk(2)*fj(1)) 569 v(5) = -sixth * (delik(1)*fi(3) + deljk(1)*fj(3) + 570 $ delik(3)*fi(1) + deljk(3)*fj(1)) 571 v(6) = -sixth * (delik(2)*fi(3) + deljk(2)*fj(3) + 572 $ delik(3)*fi(2) + deljk(3)*fj(2)) 573 574 vatom(1,i) = vatom(1,i) + v(1) 575 vatom(2,i) = vatom(2,i) + v(2) 576 vatom(3,i) = vatom(3,i) + v(3) 577 vatom(4,i) = vatom(4,i) + v(4) 578 vatom(5,i) = vatom(5,i) + v(5) 579 vatom(6,i) = vatom(6,i) + v(6) 580 vatom(1,j) = vatom(1,j) + v(1) 581 vatom(2,j) = vatom(2,j) + v(2) 582 vatom(3,j) = vatom(3,j) + v(3) 583 vatom(4,j) = vatom(4,j) + v(4) 584 vatom(5,j) = vatom(5,j) + v(5) 585 vatom(6,j) = vatom(6,j) + v(6) 586 vatom(1,k) = vatom(1,k) + v(1) 587 vatom(2,k) = vatom(2,k) + v(2) 588 vatom(3,k) = vatom(3,k) + v(3) 589 vatom(4,k) = vatom(4,k) + v(4) 590 vatom(5,k) = vatom(5,k) + v(5) 591 vatom(6,k) = vatom(6,k) + v(6) 592 endif 593 594 endif 595 endif 596c end of k loop 597 enddo 598 endif 599 100 continue 600 endif 601c end of j loop 602 enddo 603 604c else if elti=0, this is not a meam atom 605 endif 606 607 return 608 end 609