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