diff -Naur lammps-18Aug09/lib/meam/meam_data.F lammps-19Aug09/lib/meam/meam_data.F --- lammps-18Aug09/lib/meam/meam_data.F 2007-01-30 14:53:58.000000000 -0700 +++ lammps-19Aug09/lib/meam/meam_data.F 2009-08-18 11:02:09.000000000 -0600 @@ -34,6 +34,7 @@ c delr_meam = cutoff region for meam c ebound_meam = factor giving maximum boundary of sceen fcn ellipse c augt1 = flag for whether t1 coefficient should be augmented +c ialloy = flag for newer alloy formulation (as in dynamo code) c gsmooth_factor = factor determining length of G smoothing region c vind[23]D = Voight notation index maps for 2 and 3D c v2D,v3D = array of factors to apply for Voight notation @@ -65,7 +66,7 @@ real*8 Cmin_meam(maxelt,maxelt,maxelt) real*8 Cmax_meam(maxelt,maxelt,maxelt) real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt) - integer augt1 + integer augt1, ialloy real*8 gsmooth_factor integer vind2D(3,3),vind3D(3,3,3) integer v2D(6),v3D(10) diff -Naur lammps-18Aug09/lib/meam/meam_dens_final.F lammps-19Aug09/lib/meam/meam_dens_final.F --- lammps-18Aug09/lib/meam/meam_dens_final.F 2008-10-01 10:13:22.000000000 -0600 +++ lammps-19Aug09/lib/meam/meam_dens_final.F 2009-08-18 11:02:09.000000000 -0600 @@ -4,21 +4,21 @@ c int *, int *, int *, c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, double *, double *, -c double *, double *, double *, double *, int *); +c double *, double *, double *, double *, double *, int *); c c Call from pair_meam.cpp has the form: c c meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, c &eng_vdwl,eatom,ntype,type,fmap, c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], -c &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, +c &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1, c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); c subroutine meam_dens_final(nlocal, nmax, $ eflag_either, eflag_global, eflag_atom, eng_vdwl, eatom, $ ntype, type, fmap, - $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, $ Gamma, dGamma1, dGamma2, dGamma3, $ rho, rho0, rho1, rho2, rho3, fp, errorflag) @@ -29,7 +29,7 @@ integer ntype, type, fmap real*8 eng_vdwl, eatom, Arho1, Arho2 real*8 Arho2b, Arho3, Arho3b - real*8 t_ave + real*8 t_ave, tsq_ave real*8 Gamma, dGamma1, dGamma2, dGamma3 real*8 rho, rho0, rho1, rho2, rho3 real*8 fp @@ -39,6 +39,7 @@ dimension type(nmax), fmap(ntype) dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) dimension Arho3(10,nmax), Arho3b(3,nmax), t_ave(3,nmax) + dimension tsq_ave(3,nmax) dimension Gamma(nmax), dGamma1(nmax), dGamma2(nmax) dimension dGamma3(nmax), rho(nmax), rho0(nmax) dimension rho1(nmax), rho2(nmax), rho3(nmax) @@ -70,9 +71,15 @@ enddo if( rho0(i) .gt. 0.0 ) then - t_ave(1,i) = t_ave(1,i)/rho0(i) - t_ave(2,i) = t_ave(2,i)/rho0(i) - t_ave(3,i) = t_ave(3,i)/rho0(i) + if (ialloy.eq.1) then + t_ave(1,i) = t_ave(1,i)/tsq_ave(1,i) + t_ave(2,i) = t_ave(2,i)/tsq_ave(2,i) + t_ave(3,i) = t_ave(3,i)/tsq_ave(3,i) + else + t_ave(1,i) = t_ave(1,i)/rho0(i) + t_ave(2,i) = t_ave(2,i)/rho0(i) + t_ave(3,i) = t_ave(3,i)/rho0(i) + endif endif Gamma(i) = t_ave(1,i)*rho1(i) diff -Naur lammps-18Aug09/lib/meam/meam_dens_init.F lammps-19Aug09/lib/meam/meam_dens_init.F --- lammps-18Aug09/lib/meam/meam_dens_init.F 2008-10-01 10:13:22.000000000 -0600 +++ lammps-19Aug09/lib/meam/meam_dens_init.F 2009-08-18 11:02:09.000000000 -0600 @@ -3,7 +3,7 @@ c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *, c int *, int *, int *, int *, c double *, double *, double *, double *, double *, double *, -c double *, double *, double *, double *, int *); +c double *, double *, double *, double *, double *, int *); c c c Call from pair_meam.cpp has the form: @@ -12,7 +12,7 @@ c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], c rho0,&arho1[0][0],&arho2[0][0],arho2b, -c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&errorflag); +c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag); c subroutine meam_dens_init(i, nmax, @@ -20,7 +20,7 @@ $ numneigh, firstneigh, $ numneigh_full, firstneigh_full, $ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b, - $ arho3, arho3b, t_ave, errorflag) + $ arho3, arho3b, t_ave, tsq_ave, errorflag) use meam_data implicit none @@ -30,7 +30,7 @@ integer numneigh, firstneigh, numneigh_full, firstneigh_full real*8 scrfcn, dscrfcn, fcpair real*8 rho0, arho1, arho2 - real*8 arho2b, arho3, arho3b, t_ave + real*8 arho2b, arho3, arho3b, t_ave, tsq_ave integer errorflag integer j,jn @@ -40,7 +40,7 @@ dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh) dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax) dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax) - dimension t_ave(3,nmax) + dimension t_ave(3,nmax), tsq_ave(3,nmax) errorflag = 0 @@ -54,7 +54,7 @@ call calc_rho1(i, nmax, ntype, type, fmap, x, $ numneigh, firstneigh, $ scrfcn, fcpair, rho0, arho1, arho2, arho2b, - $ arho3, arho3b, t_ave) + $ arho3, arho3b, t_ave, tsq_ave) return end @@ -194,7 +194,7 @@ subroutine calc_rho1(i, nmax, ntype, type, fmap, x, $ numneigh, firstneigh, $ scrfcn, fcpair, rho0, arho1, arho2, arho2b, - $ arho3, arho3b, t_ave) + $ arho3, arho3b, t_ave, tsq_ave) use meam_data implicit none @@ -203,14 +203,14 @@ real*8 x integer numneigh, firstneigh real*8 scrfcn, fcpair, rho0, arho1, arho2 - real*8 arho2b, arho3, arho3b, t_ave + real*8 arho2b, arho3, arho3b, t_ave, tsq_ave dimension type(nmax), fmap(ntype), x(3,nmax) dimension firstneigh(numneigh) dimension scrfcn(numneigh), fcpair(numneigh) dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax) dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax) - dimension t_ave(3,nmax) + dimension t_ave(3,nmax), tsq_ave(3,nmax) integer jn,j,m,n,p,elti,eltj integer nv2,nv3 @@ -248,6 +248,14 @@ rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)*sij rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)*sij rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)*sij + if (ialloy.eq.1) then + rhoa1j = rhoa1j * t1_meam(eltj) + rhoa2j = rhoa2j * t2_meam(eltj) + rhoa3j = rhoa3j * t3_meam(eltj) + rhoa1i = rhoa1i * t1_meam(elti) + rhoa2i = rhoa2i * t2_meam(elti) + rhoa3i = rhoa3i * t3_meam(elti) + endif rho0(i) = rho0(i) + rhoa0j rho0(j) = rho0(j) + rhoa0i t_ave(1,i) = t_ave(1,i) + t1_meam(eltj)*rhoa0j @@ -256,6 +264,20 @@ t_ave(1,j) = t_ave(1,j) + t1_meam(elti)*rhoa0i t_ave(2,j) = t_ave(2,j) + t2_meam(elti)*rhoa0i t_ave(3,j) = t_ave(3,j) + t3_meam(elti)*rhoa0i + if (ialloy.eq.1) then + tsq_ave(1,i) = tsq_ave(1,i) + + $ t1_meam(eltj)*t1_meam(eltj)*rhoa0j + tsq_ave(2,i) = tsq_ave(2,i) + + $ t2_meam(eltj)*t2_meam(eltj)*rhoa0j + tsq_ave(3,i) = tsq_ave(3,i) + + $ t3_meam(eltj)*t3_meam(eltj)*rhoa0j + tsq_ave(1,j) = tsq_ave(1,j) + + $ t1_meam(elti)*t1_meam(elti)*rhoa0i + tsq_ave(2,j) = tsq_ave(2,j) + + $ t2_meam(elti)*t2_meam(elti)*rhoa0i + tsq_ave(3,j) = tsq_ave(3,j) + + $ t3_meam(elti)*t3_meam(elti)*rhoa0i + endif Arho2b(i) = Arho2b(i) + rhoa2j Arho2b(j) = Arho2b(j) + rhoa2i diff -Naur lammps-18Aug09/lib/meam/meam_force.F lammps-19Aug09/lib/meam/meam_force.F --- lammps-18Aug09/lib/meam/meam_force.F 2008-10-01 10:13:22.000000000 -0600 +++ lammps-19Aug09/lib/meam/meam_force.F 2009-08-18 11:02:09.000000000 -0600 @@ -4,7 +4,7 @@ c int *, int *, int *, int *, double *, double *, c double *, double *, double *, double *, double *, double *, c double *, double *, double *, double *, double *, double *, -c double *, double *, double *, double *, double *, int *); +c double *, double *, double *, double *, double *, double *, int *); c c Call from pair_meam.cpp has the form: c @@ -14,7 +14,7 @@ c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], -c &t_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); +c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); c subroutine meam_force(i, nmax, @@ -23,8 +23,8 @@ $ numneigh, firstneigh, numneigh_full, firstneigh_full, $ scrfcn, dscrfcn, fcpair, $ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp, - $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, f, vatom, - $ errorflag) + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f, + $ vatom, errorflag) use meam_data implicit none @@ -38,7 +38,7 @@ real*8 rho0, rho1, rho2, rho3, fp real*8 Arho1, Arho2, Arho2b real*8 Arho3, Arho3b - real*8 t_ave, f, vatom + real*8 t_ave, tsq_ave, f, vatom integer errorflag dimension eatom(nmax) @@ -50,7 +50,7 @@ dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax) dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) dimension Arho3(10,nmax), Arho3b(3,nmax) - dimension t_ave(3,nmax), f(3,nmax), vatom(6,nmax) + dimension t_ave(3,nmax), tsq_ave(3,nmax), f(3,nmax), vatom(6,nmax) integer i,j,jn,k,kn,kk,m,n,p,q integer nv2,nv3,elti,eltj,eltk,ind @@ -176,6 +176,21 @@ drhoa3j = drhoa3i endif + if (ialloy.eq.1) then + rhoa1j = rhoa1j * t1_meam(eltj) + rhoa2j = rhoa2j * t2_meam(eltj) + rhoa3j = rhoa3j * t3_meam(eltj) + rhoa1i = rhoa1i * t1_meam(elti) + rhoa2i = rhoa2i * t2_meam(elti) + rhoa3i = rhoa3i * t3_meam(elti) + drhoa1j = drhoa1j * t1_meam(eltj) + drhoa2j = drhoa2j * t2_meam(eltj) + drhoa3j = drhoa3j * t3_meam(eltj) + drhoa1i = drhoa1i * t1_meam(elti) + drhoa2i = drhoa2i * t2_meam(elti) + drhoa3i = drhoa3i * t3_meam(elti) + endif + nv2 = 1 nv3 = 1 arg1i1 = 0.d0 @@ -277,21 +292,59 @@ t2j = t_ave(2,j) t3j = t_ave(3,j) - ai = 0.d0 - if( rho0(i) .ne. 0.d0 ) then - ai = drhoa0j*sij/rho0(i) - end if - aj = 0.d0 - if( rho0(j) .ne. 0.d0 ) then - aj = drhoa0i*sij/rho0(j) - end if - - dt1dr1 = ai*(t1_meam(eltj)-t1i) - dt1dr2 = aj*(t1_meam(elti)-t1j) - dt2dr1 = ai*(t2_meam(eltj)-t2i) - dt2dr2 = aj*(t2_meam(elti)-t2j) - dt3dr1 = ai*(t3_meam(eltj)-t3i) - dt3dr2 = aj*(t3_meam(elti)-t3j) + if (ialloy.eq.1) then + + a1i = 0.d0 + a1j = 0.d0 + a2i = 0.d0 + a2j = 0.d0 + a3i = 0.d0 + a3j = 0.d0 + if ( tsq_ave(1,i) .ne. 0.d0 ) then + a1i = drhoa0j*sij/tsq_ave(1,i) + endif + if ( tsq_ave(1,j) .ne. 0.d0 ) then + a1j = drhoa0i*sij/tsq_ave(1,j) + endif + if ( tsq_ave(2,i) .ne. 0.d0 ) then + a2i = drhoa0j*sij/tsq_ave(2,i) + endif + if ( tsq_ave(2,j) .ne. 0.d0 ) then + a2j = drhoa0i*sij/tsq_ave(2,j) + endif + if ( tsq_ave(3,i) .ne. 0.d0 ) then + a3i = drhoa0j*sij/tsq_ave(3,i) + endif + if ( tsq_ave(3,j) .ne. 0.d0 ) then + a3j = drhoa0i*sij/tsq_ave(3,j) + endif + + dt1dr1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2) + dt1dr2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2) + dt2dr1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2) + dt2dr2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2) + dt3dr1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2) + dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) + + else + + ai = 0.d0 + if( rho0(i) .ne. 0.d0 ) then + ai = drhoa0j*sij/rho0(i) + end if + aj = 0.d0 + if( rho0(j) .ne. 0.d0 ) then + aj = drhoa0i*sij/rho0(j) + end if + + dt1dr1 = ai*(t1_meam(eltj)-t1i) + dt1dr2 = aj*(t1_meam(elti)-t1j) + dt2dr1 = ai*(t2_meam(eltj)-t2i) + dt2dr2 = aj*(t2_meam(elti)-t2j) + dt3dr1 = ai*(t3_meam(eltj)-t3i) + dt3dr2 = aj*(t3_meam(elti)-t3j) + + endif c Compute derivatives of total density wrt rij, sij and rij(3) call get_shpfcn(shpi,lattce_meam(elti,elti)) @@ -340,21 +393,60 @@ drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3 drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3 - ai = 0.d0 - if( rho0(i) .ne. 0.d0 ) then - ai = rhoa0j/rho0(i) - end if - aj = 0.d0 - if( rho0(j) .ne. 0.d0 ) then - aj = rhoa0i/rho0(j) - end if + if (ialloy.eq.1) then + + a1i = 0.d0 + a1j = 0.d0 + a2i = 0.d0 + a2j = 0.d0 + a3i = 0.d0 + a3j = 0.d0 + if ( tsq_ave(1,i) .ne. 0.d0 ) then + a1i = rhoa0j/tsq_ave(1,i) + endif + if ( tsq_ave(1,j) .ne. 0.d0 ) then + a1j = rhoa0i/tsq_ave(1,j) + endif + if ( tsq_ave(2,i) .ne. 0.d0 ) then + a2i = rhoa0j/tsq_ave(2,i) + endif + if ( tsq_ave(2,j) .ne. 0.d0 ) then + a2j = rhoa0i/tsq_ave(2,j) + endif + if ( tsq_ave(3,i) .ne. 0.d0 ) then + a3i = rhoa0j/tsq_ave(3,i) + endif + if ( tsq_ave(3,j) .ne. 0.d0 ) then + a3j = rhoa0i/tsq_ave(3,j) + endif + + dt1ds1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2) + dt1ds2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2) + dt2ds1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2) + dt2ds2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2) + dt3ds1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2) + dt3ds2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) + + else + + ai = 0.d0 + if( rho0(i) .ne. 0.d0 ) then + ai = rhoa0j/rho0(i) + end if + aj = 0.d0 + if( rho0(j) .ne. 0.d0 ) then + aj = rhoa0i/rho0(j) + end if + + dt1ds1 = ai*(t1_meam(eltj)-t1i) + dt1ds2 = aj*(t1_meam(elti)-t1j) + dt2ds1 = ai*(t2_meam(eltj)-t2i) + dt2ds2 = aj*(t2_meam(elti)-t2j) + dt3ds1 = ai*(t3_meam(eltj)-t3i) + dt3ds2 = aj*(t3_meam(elti)-t3j) + + endif - dt1ds1 = ai*(t1_meam(eltj)-t1i) - dt1ds2 = aj*(t1_meam(elti)-t1j) - dt2ds1 = ai*(t2_meam(eltj)-t2i) - dt2ds2 = aj*(t2_meam(elti)-t2j) - dt3ds1 = ai*(t3_meam(eltj)-t3i) - dt3ds2 = aj*(t3_meam(elti)-t3j) drhods1 = dGamma1(i)*drho0ds1 $ + dGamma2(i)* $ (dt1ds1*rho1(i)+t1i*drho1ds1 diff -Naur lammps-18Aug09/lib/meam/meam_setup_done.F lammps-19Aug09/lib/meam/meam_setup_done.F --- lammps-18Aug09/lib/meam/meam_setup_done.F 2009-03-18 12:11:39.000000000 -0600 +++ lammps-19Aug09/lib/meam/meam_setup_done.F 2009-08-18 11:02:09.000000000 -0600 @@ -158,11 +158,12 @@ integer j,a,b,nv2 real*8 astar,frac,phizbl integer n,nmax,Z1,Z2 - real*8 arat,scrn + real*8 arat,scrn,scrn2 real*8 phitmp real*8, external :: phi_meam real*8, external :: zbl + real*8, external :: compute_phi c allocate memory for array that defines the potential @@ -178,6 +179,10 @@ allocate(phirar5(nr,(neltypes*(neltypes+1))/2)) allocate(phirar6(nr,(neltypes*(neltypes+1))/2)) +c HACK for debug: compute phi_meam for Ni3Si equilibrium spacing + r = 2.47770216 + phitmp = phi_meam(r,1,2) + c loop over pairs of element types nv2 = 0 do a = 1,neltypes @@ -194,20 +199,29 @@ c if using second-nearest neighbor, solve recursive problem c (see Lee and Baskes, PRB 62(13):8564 eqn.(21)) if (nn2_meam(a,b).eq.1) then - if (a.ne.b) then - write(6,*) 'Second-neighbor currently only valid ' - write(6,*) 'if element types are the same.' +c if (a.ne.b) then +c write(6,*) 'Second-neighbor currently only valid ' +c write(6,*) 'if element types are the same.' c call error(' ') - endif +c endif call get_Zij(Z1,lattce_meam(a,b)) call get_Zij2(Z2,arat,scrn,lattce_meam(a,b), $ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) - nmax = 10 - do n = 1,nmax - phir(j,nv2) = phir(j,nv2) + - $ (-Z2*scrn/Z1)**n * phi_meam(r*arat**n,a,b) - enddo + if (lattce_meam(a,b).eq.'b1') then + phir(j,nv2) = phir(j,nv2) - + $ Z2*scrn/(2*Z1) * phi_meam(r*arat,a,a) + call get_Zij2(Z2,arat,scrn2,lattce_meam(a,b), + $ Cmin_meam(b,b,a),Cmax_meam(b,b,a)) + phir(j,nv2) = phir(j,nv2) - + $ Z2*scrn2/(2*Z1) * phi_meam(r*arat,b,b) + else + nmax = 10 + do n = 1,nmax + phir(j,nv2) = phir(j,nv2) + + $ (-Z2*scrn/Z1)**n * phi_meam(r*arat**n,a,b) + enddo + endif endif @@ -258,7 +272,7 @@ real*8 rho02,rho12,rho22,rho32 real*8 scalfac,phiaa,phibb real*8 Eu - integer Z12 + integer Z12, errorflag character*3 latta,lattb real*8, external :: erose @@ -285,8 +299,9 @@ else c average weighting factors for the reference structure, eqn. I.8 call get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, - $ t1_meam(a),t2_meam(a),t3_meam(a), - $ t1_meam(b),t2_meam(b),t3_meam(b),lattce_meam(a,b)) + $ t1_meam(a),t2_meam(a),t3_meam(a), + $ t1_meam(b),t2_meam(b),t3_meam(b), + $ r,a,b,lattce_meam(a,b)) endif c for c11b structure, calculate background electron densities @@ -316,6 +331,12 @@ c -- note: The shape factors for the individual elements are used, since c this function should cancel the F(rho) terms when the atoms are c in the reference structure. +c -- GJW (6/12/09): I suspect there may be a problem here... since we should be +c using the pure element reference structure, we should also +c be using the element "t" values (not the averages compute +c above). It doesn't matter for reference structures with +c s(1)=s(2)=s(3)=0, but might cause diffs otherwise. For now +c what's here is consistent with what's done in meam_dens_final. Z1 = Z_meam(a) Z2 = Z_meam(b) if (ibar_meam(a).le.0) then @@ -323,14 +344,14 @@ else call get_shpfcn(s1,lattce_meam(a,a)) Gam1 = (s1(1)*t11av+s1(2)*t21av+s1(3)*t31av)/(Z1*Z1) - call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1) + call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag) endif if (ibar_meam(b).le.0) then G2 = 1.d0 else call get_shpfcn(s2,lattce_meam(b,b)) Gam2 = (s2(1)*t12av+s2(2)*t22av+s2(3)*t32av)/(Z2*Z2) - call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2) + call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag) endif rho0_1 = rho0_meam(a)*Z1*G1 rho0_2 = rho0_meam(b)*Z2*G2 @@ -347,8 +368,8 @@ Gam1 = Gam1/(rho01*rho01) Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32) Gam2 = Gam2/(rho02*rho02) - call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1) - call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2) + call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag) + call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag) rhobar1 = rho01/rho0_1*G1 rhobar2 = rho02/rho0_2*G2 endif @@ -370,17 +391,20 @@ c calculate the pair energy if (lattce_meam(a,b).eq.'c11') then - latta = lattce_meam(a,a) - if (latta.eq.'dia') then - phiaa = phi_meam(r,a,a) - phi_m = (3*Eu - F2 - 2*F1 - 5*phiaa)/Z12 - else - phibb = phi_meam(r,b,b) - phi_m = (3*Eu - F1 - 2*F2 - 5*phibb)/Z12 - endif -c phi_m = 0.d0 + latta = lattce_meam(a,a) + if (latta.eq.'dia') then + phiaa = phi_meam(r,a,a) + phi_m = (3*Eu - F2 - 2*F1 - 5*phiaa)/Z12 + else + phibb = phi_meam(r,b,b) + phi_m = (3*Eu - F1 - 2*F2 - 5*phibb)/Z12 + endif +c phi_m = 0.d0 + else if (lattce_meam(a,b).eq.'l12') then + phiaa = phi_meam(r,a,a) + phi_m = Eu/3. - F1/4. - F2/12. - phiaa else -c +c c potential is computed from Rose function and embedding energy phi_m = (2*Eu - F1 - F2)/Z12 c @@ -426,11 +450,15 @@ c------------------------------------------------------------------------------c c Average weighting factors for the reference structure subroutine get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, - $ t11,t21,t31,t12,t22,t32,latt) + $ t11,t21,t31,t12,t22,t32, + $ r,a,b,latt) + use meam_data implicit none real*8 t11av,t21av,t31av,t12av,t22av,t32av - real*8 t11,t21,t31,t12,t22,t32 + real*8 t11,t21,t31,t12,t22,t32,r + integer a,b character*3 latt + real*8 rhoa01,rhoa02,a1,a2,rho01,rho02 if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'dia' $ .or.latt.eq.'hcp'.or.latt.eq.'b1' $ .or.latt.eq.'dim') then @@ -442,7 +470,21 @@ t22av = t21 t32av = t31 else -c call error('Lattice not defined in get_tavref.') + a1 = r/re_meam(a,a) - 1.d0 + a2 = r/re_meam(b,b) - 1.d0 + rhoa01 = rho0_meam(a)*exp(-beta0_meam(a)*a1) + rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2) + if (latt.eq.'l12') then + rho01 = 8*rhoa01 + 4*rhoa02 + t11av = (8*t11*rhoa01 + 4*t12*rhoa02)/rho01 + t12av = t11 + t21av = (8*t21*rhoa01 + 4*t22*rhoa02)/rho01 + t22av = t21 + t31av = (8*t31*rhoa01 + 4*t32*rhoa02)/rho01 + t32av = t31 + else +c call error('Lattice not defined in get_tavref.') + endif endif return end @@ -466,6 +508,8 @@ Zij = 1 else if (latt.eq.'c11') then Zij = 10 + else if (latt.eq.'l12') then + Zij = 12 else c call error('Lattice not defined in get_Zij.') endif @@ -504,9 +548,9 @@ a = sqrt(2.d0) numscr = 4 else if (latt.eq.'b1') then - Zij2 = 6 + Zij2 = 12 a = sqrt(2.d0) - numscr = 4 + numscr = 2 else c call error('Lattice not defined in get_Zij2.') endif @@ -535,7 +579,7 @@ integer a,b integer Zij1nn,Zij2nn real*8 rhoa01nn,rhoa02nn - real*8 arat,scrn + real*8 arat,scrn,denom a1 = r/re_meam(a,a) - 1.d0 a2 = r/re_meam(b,b) - 1.d0 @@ -567,17 +611,6 @@ rho01 = 8.d0*rhoa02 rho02 = 8.d0*rhoa01 else if (lat.eq.'b1') then -c Warning: this assumes that atoms of element 1 are not completely -c screened from each other, but that atoms of element 2 are. So -c Cmin_meam(1,1,2) can be less than 1.0, but Cmin_meam(2,2,1) cannot. -c This should be made more general in the future. -c Cmin = Cmin_meam(1,1,2) -c Cmax = Cmax_meam(1,1,2) -c C = (1.0-Cmin)/(Cmax-Cmin) -c call fcut(C,fc) -c a1 = sqrt(2.0)*r/re_meam(1,1) - 1.d0 -c rho01 = 6*rhoa02 + 12*fc*fc*rho0_meam(1)*exp(-beta0_meam(1)*a1) -c rho02 = 6*rhoa01 rho01 = 6*rhoa02 rho02 = 6*rhoa01 else if (lat.eq.'dia') then @@ -609,6 +642,18 @@ rho22 = rhoa22 rho31 = rhoa31 rho32 = rhoa32 + else if (lat.eq.'l12') then + rho01 = 8*rhoa01 + 4*rhoa02 + rho02 = 12*rhoa01 + if (ialloy.eq.0) then + rho21 = 8./3.*(rhoa21-rhoa22)*(rhoa21-rhoa22) + else + rho21 = 8./3.*(rhoa21*t2_meam(a)-rhoa22*t2_meam(b))**2 + denom = 8*rhoa01*t2_meam(a)**2 + 4*rhoa02*t2_meam(b)**2 + if (denom.gt.0.) then + rho21 = rho21/denom * rho01 + endif + endif else c call error('Lattice not defined in get_densref.') endif @@ -617,9 +662,9 @@ c Assume that second neighbor is of same type, first neighbor c may be of different type - if (lat.ne.'bcc') then +c if (lat.ne.'bcc') then c call error('Second-neighbor not defined for this lattice') - endif +c endif call get_Zij2(Zij2nn,arat,scrn,lat, $ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) @@ -630,8 +675,12 @@ rhoa01nn = rho0_meam(a)*exp(-beta0_meam(a)*a1) rhoa02nn = rho0_meam(b)*exp(-beta0_meam(b)*a2) - rho01 = rho01 + Zij2nn*scrn*rhoa02nn - rho02 = rho02 + Zij2nn*scrn*rhoa01nn + rho01 = rho01 + Zij2nn*scrn*rhoa01nn + +c Assume Zij2nn and arat don't depend on order, but scrn might + call get_Zij2(Zij2nn,arat,scrn,lat, + $ Cmin_meam(b,b,a),Cmax_meam(b,b,a)) + rho02 = rho02 + Zij2nn*scrn*rhoa02nn endif @@ -731,3 +780,25 @@ enddo end + +c--------------------------------------------------------------------- +c Compute Rose energy function, I.16 +c + real*8 function compute_phi(rij, elti, eltj) + use meam_data + implicit none + + real*8 rij, pp + integer elti, eltj, ind, kk + + ind = eltind(elti, eltj) + pp = rij*rdrar + 1.0D0 + kk = pp + kk = min(kk,nrar-1) + pp = pp - kk + pp = min(pp,1.0D0) + compute_phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp + $ + phirar1(kk,ind))*pp + phirar(kk,ind) + + return + end diff -Naur lammps-18Aug09/lib/meam/meam_setup_global.F lammps-19Aug09/lib/meam/meam_setup_global.F --- lammps-18Aug09/lib/meam/meam_setup_global.F 2007-01-30 14:53:58.000000000 -0700 +++ lammps-19Aug09/lib/meam/meam_setup_global.F 2009-08-18 11:02:09.000000000 -0600 @@ -98,6 +98,7 @@ nn2_meam(:,:) = 0 gsmooth_factor = 99.0 augt1 = 1 + ialloy = 0 return end diff -Naur lammps-18Aug09/lib/meam/meam_setup_param.F lammps-19Aug09/lib/meam/meam_setup_param.F --- lammps-18Aug09/lib/meam/meam_setup_param.F 2007-01-30 14:53:58.000000000 -0700 +++ lammps-19Aug09/lib/meam/meam_setup_param.F 2009-08-18 11:02:09.000000000 -0600 @@ -27,6 +27,7 @@ c 12 = augt1 c 13 = gsmooth_factor c 14 = re_meam +c 15 = ialloy subroutine meam_setup_param(which, value, nindex, $ index, errorflag) @@ -71,6 +72,8 @@ lattce_meam(index(1),index(2)) = 'b1' else if (value.eq.6) then lattce_meam(index(1),index(2)) = 'c11' + else if (value.eq.7) then + lattce_meam(index(1),index(2)) = 'l12' endif c 5 = attrac_meam @@ -113,6 +116,10 @@ else if (which.eq.14) then re_meam(index(1),index(2)) = value +c 15 = ialloy + else if (which.eq.15) then + ialloy = value + else errorflag = 1 endif diff -Naur lammps-18Aug09/src/MEAM/pair_meam.cpp lammps-19Aug09/src/MEAM/pair_meam.cpp --- lammps-18Aug09/src/MEAM/pair_meam.cpp 2008-02-19 09:23:58.000000000 -0700 +++ lammps-19Aug09/src/MEAM/pair_meam.cpp 2009-08-18 11:01:49.000000000 -0600 @@ -1,901 +1,919 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Greg Wagner (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_meam.h" -#include "atom.h" -#include "force.h" -#include "comm.h" -#include "memory.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - -#define MAXLINE 1024 - -enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11}; -int nkeywords = 15; -char *keywords[] = {"Ec","alpha","rho0","delta","lattce", - "attrac","repuls","nn2","Cmin","Cmax","rc","delr", - "augt1","gsmooth_factor","re"}; - -/* ---------------------------------------------------------------------- */ - -PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 0; - one_coeff = 1; - - nmax = 0; - rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL; - gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL; - arho1 = arho2 = arho3 = arho3b = t_ave = NULL; - - maxneigh = 0; - scrfcn = dscrfcn = fcpair = NULL; - - nelements = 0; - elements = NULL; - mass = NULL; - - // set comm size needed by this Pair - - comm_forward = 35; - comm_reverse = 27; -} - -/* ---------------------------------------------------------------------- - free all arrays - check if allocated, since class can be destructed when incomplete -------------------------------------------------------------------------- */ - -PairMEAM::~PairMEAM() -{ - meam_cleanup_(); - - memory->sfree(rho); - memory->sfree(rho0); - memory->sfree(rho1); - memory->sfree(rho2); - memory->sfree(rho3); - memory->sfree(frhop); - memory->sfree(gamma); - memory->sfree(dgamma1); - memory->sfree(dgamma2); - memory->sfree(dgamma3); - memory->sfree(arho2b); - - memory->destroy_2d_double_array(arho1); - memory->destroy_2d_double_array(arho2); - memory->destroy_2d_double_array(arho3); - memory->destroy_2d_double_array(arho3b); - memory->destroy_2d_double_array(t_ave); - - memory->sfree(scrfcn); - memory->sfree(dscrfcn); - memory->sfree(fcpair); - - for (int i = 0; i < nelements; i++) delete [] elements[i]; - delete [] elements; - delete [] mass; - - if (allocated) { - memory->destroy_2d_int_array(setflag); - memory->destroy_2d_double_array(cutsq); - delete [] map; - delete [] fmap; - } -} - -/* ---------------------------------------------------------------------- */ - -void PairMEAM::compute(int eflag, int vflag) -{ - int i,j,ii,n,inum_half,itype,jtype,errorflag; - double evdwl; - int *ilist_half,*jlist_half,*numneigh_half,**firstneigh_half; - int *numneigh_full,**firstneigh_full; - - evdwl = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int newton_pair = force->newton_pair; - - // grow local arrays if necessary - - if (atom->nmax > nmax) { - memory->sfree(rho); - memory->sfree(rho0); - memory->sfree(rho1); - memory->sfree(rho2); - memory->sfree(rho3); - memory->sfree(frhop); - memory->sfree(gamma); - memory->sfree(dgamma1); - memory->sfree(dgamma2); - memory->sfree(dgamma3); - memory->sfree(arho2b); - memory->destroy_2d_double_array(arho1); - memory->destroy_2d_double_array(arho2); - memory->destroy_2d_double_array(arho3); - memory->destroy_2d_double_array(arho3b); - memory->destroy_2d_double_array(t_ave); - - nmax = atom->nmax; - - rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho"); - rho0 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho0"); - rho1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho1"); - rho2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho2"); - rho3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho3"); - frhop = (double *) memory->smalloc(nmax*sizeof(double),"pair:frhop"); - gamma = (double *) memory->smalloc(nmax*sizeof(double),"pair:gamma"); - dgamma1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma1"); - dgamma2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma2"); - dgamma3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma3"); - arho2b = (double *) memory->smalloc(nmax*sizeof(double),"pair:arho2b"); - arho1 = memory->create_2d_double_array(nmax,3,"pair:arho1"); - arho2 = memory->create_2d_double_array(nmax,6,"pair:arho2"); - arho3 = memory->create_2d_double_array(nmax,10,"pair:arho3"); - arho3b = memory->create_2d_double_array(nmax,3,"pair:arho3b"); - t_ave = memory->create_2d_double_array(nmax,3,"pair:t_ave"); - } - - // neighbor list info - - inum_half = listhalf->inum; - ilist_half = listhalf->ilist; - numneigh_half = listhalf->numneigh; - firstneigh_half = listhalf->firstneigh; - numneigh_full = listfull->numneigh; - firstneigh_full = listfull->firstneigh; - - // check size of scrfcn based on half neighbor list - - int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - - n = 0; - for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]]; - - if (n > maxneigh) { - memory->sfree(scrfcn); - memory->sfree(dscrfcn); - memory->sfree(fcpair); - maxneigh = n; - scrfcn = - (double *) memory->smalloc(maxneigh*sizeof(double),"pair:scrfcn"); - dscrfcn = - (double *) memory->smalloc(maxneigh*sizeof(double),"pair:dscrfcn"); - fcpair = - (double *) memory->smalloc(maxneigh*sizeof(double),"pair:fcpair"); - } - - // zero out local arrays - - for (i = 0; i < nall; i++) { - rho0[i] = 0.0; - arho2b[i] = 0.0; - arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0; - for (j = 0; j < 6; j++) arho2[i][j] = 0.0; - for (j = 0; j < 10; j++) arho3[i][j] = 0.0; - arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0; - t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0; - } - - double **x = atom->x; - double **f = atom->f; - int *type = atom->type; - int ntype = atom->ntypes; - - // change neighbor list indices to Fortran indexing - - neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half); - neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full); - - // 3 stages of MEAM calculation - // loop over my atoms followed by communication - - int ifort; - int offset = 0; - errorflag = 0; - - for (ii = 0; ii < inum_half; ii++) { - i = ilist_half[ii]; - ifort = i+1; - meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0], - &numneigh_half[i],firstneigh_half[i], - &numneigh_full[i],firstneigh_full[i], - &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], - rho0,&arho1[0][0],&arho2[0][0],arho2b, - &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&errorflag); - if (errorflag) { - char str[128]; - sprintf(str,"MEAM library error %d",errorflag); - error->one(str); - } - offset += numneigh_half[i]; - } - - comm->reverse_comm_pair(this); - - meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, - &eng_vdwl,eatom,&ntype,type,fmap, - &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], - &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, - dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); - if (errorflag) { - char str[128]; - sprintf(str,"MEAM library error %d",errorflag); - error->one(str); - } - - comm->comm_pair(this); - - offset = 0; - - // vptr is first value in vatom if it will be used by meam_force() - // else vatom may not exist, so pass dummy ptr - - double *vptr; - if (vflag_atom) vptr = &vatom[0][0]; - else vptr = &cutmax; - - for (ii = 0; ii < inum_half; ii++) { - i = ilist_half[ii]; - ifort = i+1; - meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, - &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], - &numneigh_half[i],firstneigh_half[i], - &numneigh_full[i],firstneigh_full[i], - &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], - dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, - &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], - &t_ave[0][0],&f[0][0],vptr,&errorflag); - if (errorflag) { - char str[128]; - sprintf(str,"MEAM library error %d",errorflag); - error->one(str); - } - offset += numneigh_half[i]; - } - - // change neighbor list indices back to C indexing - - neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); - neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); - - if (vflag_fdotr) virial_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairMEAM::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); - cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); - - map = new int[n+1]; - fmap = new int[n]; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairMEAM::settings(int narg, char **arg) -{ - if (narg != 0) error->all("Illegal pair_style command"); -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairMEAM::coeff(int narg, char **arg) -{ - int i,j,m,n; - - if (!allocated) allocate(); - - if (narg < 6) error->all("Incorrect args for pair coefficients"); - - // insure I,J args are * * - - if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) - error->all("Incorrect args for pair coefficients"); - - // read MEAM element names between 2 filenames - // nelements = # of MEAM elements - // elements = list of unique element names - - if (nelements) { - for (i = 0; i < nelements; i++) delete [] elements[i]; - delete [] elements; - delete [] mass; - } - nelements = narg - 4 - atom->ntypes; - if (nelements < 1) error->all("Incorrect args for pair coefficients"); - elements = new char*[nelements]; - mass = new double[nelements]; - - for (i = 0; i < nelements; i++) { - n = strlen(arg[i+3]) + 1; - elements[i] = new char[n]; - strcpy(elements[i],arg[i+3]); - } - - // read MEAM library and parameter files - // pass all parameters to MEAM package - // tell MEAM package that setup is done - - read_files(arg[2],arg[2+nelements+1]); - meam_setup_done_(&cutmax); - - // read args that map atom types to MEAM elements - // map[i] = which element the Ith atom type is, -1 if not mapped - - for (i = 4 + nelements; i < narg; i++) { - m = i - (4+nelements) + 1; - for (j = 0; j < nelements; j++) - if (strcmp(arg[i],elements[j]) == 0) break; - if (j < nelements) map[m] = j; - else if (strcmp(arg[i],"NULL") == 0) map[m] = -1; - else error->all("Incorrect args for pair coefficients"); - } - - // clear setflag since coeff() called once with I,J = * * - - n = atom->ntypes; - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - // set setflag i,j for type pairs where both are mapped to elements - // set mass for i,i in atom class - - int count = 0; - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - if (map[i] >= 0 && map[j] >= 0) { - setflag[i][j] = 1; - if (i == j) atom->set_mass(i,mass[map[i]]); - count++; - } - - if (count == 0) error->all("Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairMEAM::init_style() -{ - if (force->newton_pair == 0) - error->all("Pair style MEAM requires newton pair on"); - - // need full and half neighbor list - - int irequest_full = neighbor->request(this); - neighbor->requests[irequest_full]->id = 1; - neighbor->requests[irequest_full]->half = 0; - neighbor->requests[irequest_full]->full = 1; - int irequest_half = neighbor->request(this); - neighbor->requests[irequest_half]->id = 2; - neighbor->requests[irequest_half]->half = 0; - neighbor->requests[irequest_half]->half_from_full = 1; - neighbor->requests[irequest_half]->otherlist = irequest_full; - - // setup Fortran-style mapping array needed by MEAM package - // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index - // if type I is not a MEAM atom, fmap stores a 0 - - for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1; -} - -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - half or full -------------------------------------------------------------------------- */ - -void PairMEAM::init_list(int id, NeighList *ptr) -{ - if (id == 1) listfull = ptr; - else if (id == 2) listhalf = ptr; -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairMEAM::init_one(int i, int j) -{ - return cutmax; -} - -/* ---------------------------------------------------------------------- */ - -void PairMEAM::read_files(char *globalfile, char *userfile) -{ - // open global meamf file on proc 0 - - FILE *fp; - if (comm->me == 0) { - fp = fopen(globalfile,"r"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open MEAM potential file %s",globalfile); - error->one(str); - } - } - - // allocate parameter arrays - - int params_per_line = 19; - - int *lat = new int[nelements]; - int *ielement = new int[nelements]; - int *ibar = new int[nelements]; - double *z = new double[nelements]; - double *atwt = new double[nelements]; - double *alpha = new double[nelements]; - double *b0 = new double[nelements]; - double *b1 = new double[nelements]; - double *b2 = new double[nelements]; - double *b3 = new double[nelements]; - double *alat = new double[nelements]; - double *esub = new double[nelements]; - double *asub = new double[nelements]; - double *t0 = new double[nelements]; - double *t1 = new double[nelements]; - double *t2 = new double[nelements]; - double *t3 = new double[nelements]; - double *rozero = new double[nelements]; - - bool *found = new bool[nelements]; - for (int i = 0; i < nelements; i++) found[i] = false; - - // read each set of params from global MEAM file - // one set of params can span multiple lines - // store params if element name is in element list - // if element name appears multiple times, only store 1st entry - - int i,n,nwords; - char **words = new char*[params_per_line+1]; - char line[MAXLINE],*ptr; - int eof = 0; - - int nset = 0; - while (1) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fp); - if (ptr == NULL) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if (ptr = strchr(line,'#')) *ptr = '\0'; - nwords = atom->count_words(line); - if (nwords == 0) continue; - - // concatenate additional lines until have params_per_line words - - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n],MAXLINE-n,fp); - if (ptr == NULL) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - if (ptr = strchr(line,'#')) *ptr = '\0'; - nwords = atom->count_words(line); - } - - if (nwords != params_per_line) - error->all("Incorrect format in MEAM potential file"); - - // words = ptrs to all words in line - // strip single and double quotes from words - - nwords = 0; - words[nwords++] = strtok(line,"' \t\n\r\f"); - while (words[nwords++] = strtok(NULL,"' \t\n\r\f")) continue; - - // skip if element name isn't in element list - - for (i = 0; i < nelements; i++) - if (strcmp(words[0],elements[i]) == 0) break; - if (i == nelements) continue; - - // skip if element already appeared - - if (found[i] == true) continue; - found[i] = true; - - // map lat string to an integer - - if (strcmp(words[1],"fcc") == 0) lat[i] = FCC; - else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC; - else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP; - else if (strcmp(words[1],"dim") == 0) lat[i] = DIM; - else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND; - else error->all("Unrecognized lattice type in MEAM file 1"); - - // store parameters - - z[i] = atof(words[2]); - ielement[i] = atoi(words[3]); - atwt[i] = atof(words[4]); - alpha[i] = atof(words[5]); - b0[i] = atof(words[6]); - b1[i] = atof(words[7]); - b2[i] = atof(words[8]); - b3[i] = atof(words[9]); - alat[i] = atof(words[10]); - esub[i] = atof(words[11]); - asub[i] = atof(words[12]); - t0[i] = atof(words[13]); - t1[i] = atof(words[14]); - t2[i] = atof(words[15]); - t3[i] = atof(words[16]); - rozero[i] = atof(words[17]); - ibar[i] = atoi(words[18]); - - nset++; - } - - // error if didn't find all elements in file - - if (nset != nelements) - error->all("Did not find all elements in MEAM library file"); - - // pass element parameters to MEAM package - - meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, - alat,esub,asub,t0,t1,t2,t3,rozero,ibar); - - // set element masses - - for (i = 0; i < nelements; i++) mass[i] = atwt[i]; - - // clean-up memory - - delete [] words; - - delete [] lat; - delete [] ielement; - delete [] ibar; - delete [] z; - delete [] atwt; - delete [] alpha; - delete [] b0; - delete [] b1; - delete [] b2; - delete [] b3; - delete [] alat; - delete [] esub; - delete [] asub; - delete [] t0; - delete [] t1; - delete [] t2; - delete [] t3; - delete [] rozero; - delete [] found; - - // done if user param file is NULL - - if (strcmp(userfile,"NULL") == 0) return; - - // open user param file on proc 0 - - if (comm->me == 0) { - fp = fopen(userfile,"r"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open MEAM potential file %s",userfile); - error->one(str); - } - } - - // read settings - // pass them one at a time to MEAM package - // match strings to list of corresponding ints - - int which; - double value; - int nindex,index[3]; - int maxparams = 6; - char **params = new char*[maxparams]; - int nparams; - - eof = 0; - while (1) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fp); - if (ptr == NULL) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; - } - MPI_Bcast(&eof,1,MPI_INT,0,world); - if (eof) break; - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - - // strip comment, skip line if blank - - if (ptr = strchr(line,'#')) *ptr = '\0'; - nparams = atom->count_words(line); - if (nparams == 0) continue; - - // words = ptrs to all words in line - - nparams = 0; - params[nparams++] = strtok(line,"=(), '\t\n\r\f"); - while (nparams < maxparams && - (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f"))) - continue; - nparams--; - - for (which = 0; which < nkeywords; which++) - if (strcmp(params[0],keywords[which]) == 0) break; - if (which == nkeywords) { - char str[128]; - sprintf(str,"Keyword %s in MEAM parameter file not recognized", - params[0]); - error->all(str); - } - nindex = nparams - 2; - for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]); - - // map lattce_meam value to an integer - - if (which == 4) { - if (strcmp(params[nparams-1],"fcc") == 0) value = FCC; - else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC; - else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP; - else if (strcmp(params[nparams-1],"dim") == 0) value = DIM; - else if (strcmp(params[nparams-1],"dia") == 0) value = DIAMOND; - else if (strcmp(params[nparams-1],"b1") == 0) value = B1; - else if (strcmp(params[nparams-1],"c11") == 0) value = C11; - else error->all("Unrecognized lattice type in MEAM file 2"); - } - else value = atof(params[nparams-1]); - - // pass single setting to MEAM package - - int errorflag = 0; - meam_setup_param_(&which,&value,&nindex,index,&errorflag); - if (errorflag) { - char str[128]; - sprintf(str,"MEAM library error %d",errorflag); - error->all(str); - } - } - - delete [] params; -} - -/* ---------------------------------------------------------------------- */ - -int PairMEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) -{ - int i,j,k,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = rho0[j]; - buf[m++] = rho1[j]; - buf[m++] = rho2[j]; - buf[m++] = rho3[j]; - buf[m++] = frhop[j]; - buf[m++] = gamma[j]; - buf[m++] = dgamma1[j]; - buf[m++] = dgamma2[j]; - buf[m++] = dgamma3[j]; - buf[m++] = arho2b[j]; - buf[m++] = arho1[j][0]; - buf[m++] = arho1[j][1]; - buf[m++] = arho1[j][2]; - buf[m++] = arho2[j][0]; - buf[m++] = arho2[j][1]; - buf[m++] = arho2[j][2]; - buf[m++] = arho2[j][3]; - buf[m++] = arho2[j][4]; - buf[m++] = arho2[j][5]; - for (k = 0; k < 10; k++) buf[m++] = arho3[j][k]; - buf[m++] = arho3b[j][0]; - buf[m++] = arho3b[j][1]; - buf[m++] = arho3b[j][2]; - buf[m++] = t_ave[j][0]; - buf[m++] = t_ave[j][1]; - buf[m++] = t_ave[j][2]; - } - return 35; -} - -/* ---------------------------------------------------------------------- */ - -void PairMEAM::unpack_comm(int n, int first, double *buf) -{ - int i,k,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - rho0[i] = buf[m++]; - rho1[i] = buf[m++]; - rho2[i] = buf[m++]; - rho3[i] = buf[m++]; - frhop[i] = buf[m++]; - gamma[i] = buf[m++]; - dgamma1[i] = buf[m++]; - dgamma2[i] = buf[m++]; - dgamma3[i] = buf[m++]; - arho2b[i] = buf[m++]; - arho1[i][0] = buf[m++]; - arho1[i][1] = buf[m++]; - arho1[i][2] = buf[m++]; - arho2[i][0] = buf[m++]; - arho2[i][1] = buf[m++]; - arho2[i][2] = buf[m++]; - arho2[i][3] = buf[m++]; - arho2[i][4] = buf[m++]; - arho2[i][5] = buf[m++]; - for (k = 0; k < 10; k++) arho3[i][k] = buf[m++]; - arho3b[i][0] = buf[m++]; - arho3b[i][1] = buf[m++]; - arho3b[i][2] = buf[m++]; - t_ave[i][0] = buf[m++]; - t_ave[i][1] = buf[m++]; - t_ave[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int PairMEAM::pack_reverse_comm(int n, int first, double *buf) -{ - int i,k,m,last,size; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = rho0[i]; - buf[m++] = arho2b[i]; - buf[m++] = arho1[i][0]; - buf[m++] = arho1[i][1]; - buf[m++] = arho1[i][2]; - buf[m++] = arho2[i][0]; - buf[m++] = arho2[i][1]; - buf[m++] = arho2[i][2]; - buf[m++] = arho2[i][3]; - buf[m++] = arho2[i][4]; - buf[m++] = arho2[i][5]; - for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; - buf[m++] = arho3b[i][0]; - buf[m++] = arho3b[i][1]; - buf[m++] = arho3b[i][2]; - buf[m++] = t_ave[i][0]; - buf[m++] = t_ave[i][1]; - buf[m++] = t_ave[i][2]; - } - - return 27; -} - -/* ---------------------------------------------------------------------- */ - -void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf) -{ - int i,j,k,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - rho0[j] += buf[m++]; - arho2b[j] += buf[m++]; - arho1[j][0] += buf[m++]; - arho1[j][1] += buf[m++]; - arho1[j][2] += buf[m++]; - arho2[j][0] += buf[m++]; - arho2[j][1] += buf[m++]; - arho2[j][2] += buf[m++]; - arho2[j][3] += buf[m++]; - arho2[j][4] += buf[m++]; - arho2[j][5] += buf[m++]; - for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; - arho3b[j][0] += buf[m++]; - arho3b[j][1] += buf[m++]; - arho3b[j][2] += buf[m++]; - t_ave[j][0] += buf[m++]; - t_ave[j][1] += buf[m++]; - t_ave[j][2] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ - -double PairMEAM::memory_usage() -{ - double bytes = 11 * nmax * sizeof(double); - bytes += (3 + 6 + 10 + 3 + 3) * nmax * sizeof(double); - bytes += 3 * maxneigh * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - toggle neighbor list indices between zero- and one-based values - needed for access by MEAM Fortran library -------------------------------------------------------------------------- */ - -void PairMEAM::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh) -{ - int i,j,ii,jnum; - int *jlist; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - for (j = 0; j < jnum; j++) jlist[j]--; - } -} - -void PairMEAM::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh) -{ - int i,j,ii,jnum; - int *jlist; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - for (j = 0; j < jnum; j++) jlist[j]++; - } -} +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Greg Wagner (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_meam.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +#define MAXLINE 1024 + +enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11,L12}; +int nkeywords = 16; +char *keywords[] = {"Ec","alpha","rho0","delta","lattce", + "attrac","repuls","nn2","Cmin","Cmax","rc","delr", + "augt1","gsmooth_factor","re","ialloy"}; + +/* ---------------------------------------------------------------------- */ + +PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + one_coeff = 1; + + nmax = 0; + rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL; + gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL; + arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL; + + maxneigh = 0; + scrfcn = dscrfcn = fcpair = NULL; + + nelements = 0; + elements = NULL; + mass = NULL; + + // set comm size needed by this Pair + + comm_forward = 35; + comm_reverse = 27; +} + +/* ---------------------------------------------------------------------- + free all arrays + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairMEAM::~PairMEAM() +{ + meam_cleanup_(); + + memory->sfree(rho); + memory->sfree(rho0); + memory->sfree(rho1); + memory->sfree(rho2); + memory->sfree(rho3); + memory->sfree(frhop); + memory->sfree(gamma); + memory->sfree(dgamma1); + memory->sfree(dgamma2); + memory->sfree(dgamma3); + memory->sfree(arho2b); + + memory->destroy_2d_double_array(arho1); + memory->destroy_2d_double_array(arho2); + memory->destroy_2d_double_array(arho3); + memory->destroy_2d_double_array(arho3b); + memory->destroy_2d_double_array(t_ave); + memory->destroy_2d_double_array(tsq_ave); + + memory->sfree(scrfcn); + memory->sfree(dscrfcn); + memory->sfree(fcpair); + + for (int i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + delete [] mass; + + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + delete [] map; + delete [] fmap; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairMEAM::compute(int eflag, int vflag) +{ + int i,j,ii,n,inum_half,itype,jtype,errorflag; + double evdwl; + int *ilist_half,*jlist_half,*numneigh_half,**firstneigh_half; + int *numneigh_full,**firstneigh_full; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + int newton_pair = force->newton_pair; + + // grow local arrays if necessary + + if (atom->nmax > nmax) { + memory->sfree(rho); + memory->sfree(rho0); + memory->sfree(rho1); + memory->sfree(rho2); + memory->sfree(rho3); + memory->sfree(frhop); + memory->sfree(gamma); + memory->sfree(dgamma1); + memory->sfree(dgamma2); + memory->sfree(dgamma3); + memory->sfree(arho2b); + memory->destroy_2d_double_array(arho1); + memory->destroy_2d_double_array(arho2); + memory->destroy_2d_double_array(arho3); + memory->destroy_2d_double_array(arho3b); + memory->destroy_2d_double_array(t_ave); + memory->destroy_2d_double_array(tsq_ave); + + nmax = atom->nmax; + + rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho"); + rho0 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho0"); + rho1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho1"); + rho2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho2"); + rho3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho3"); + frhop = (double *) memory->smalloc(nmax*sizeof(double),"pair:frhop"); + gamma = (double *) memory->smalloc(nmax*sizeof(double),"pair:gamma"); + dgamma1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma1"); + dgamma2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma2"); + dgamma3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma3"); + arho2b = (double *) memory->smalloc(nmax*sizeof(double),"pair:arho2b"); + arho1 = memory->create_2d_double_array(nmax,3,"pair:arho1"); + arho2 = memory->create_2d_double_array(nmax,6,"pair:arho2"); + arho3 = memory->create_2d_double_array(nmax,10,"pair:arho3"); + arho3b = memory->create_2d_double_array(nmax,3,"pair:arho3b"); + t_ave = memory->create_2d_double_array(nmax,3,"pair:t_ave"); + tsq_ave = memory->create_2d_double_array(nmax,3,"pair:tsq_ave"); + } + + // neighbor list info + + inum_half = listhalf->inum; + ilist_half = listhalf->ilist; + numneigh_half = listhalf->numneigh; + firstneigh_half = listhalf->firstneigh; + numneigh_full = listfull->numneigh; + firstneigh_full = listfull->firstneigh; + + // check size of scrfcn based on half neighbor list + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + n = 0; + for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]]; + + if (n > maxneigh) { + memory->sfree(scrfcn); + memory->sfree(dscrfcn); + memory->sfree(fcpair); + maxneigh = n; + scrfcn = + (double *) memory->smalloc(maxneigh*sizeof(double),"pair:scrfcn"); + dscrfcn = + (double *) memory->smalloc(maxneigh*sizeof(double),"pair:dscrfcn"); + fcpair = + (double *) memory->smalloc(maxneigh*sizeof(double),"pair:fcpair"); + } + + // zero out local arrays + + for (i = 0; i < nall; i++) { + rho0[i] = 0.0; + arho2b[i] = 0.0; + arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0; + for (j = 0; j < 6; j++) arho2[i][j] = 0.0; + for (j = 0; j < 10; j++) arho3[i][j] = 0.0; + arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0; + t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0; + tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0; + } + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int ntype = atom->ntypes; + + // change neighbor list indices to Fortran indexing + + neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half); + neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full); + + // 3 stages of MEAM calculation + // loop over my atoms followed by communication + + int ifort; + int offset = 0; + errorflag = 0; + + for (ii = 0; ii < inum_half; ii++) { + i = ilist_half[ii]; + ifort = i+1; + meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0], + &numneigh_half[i],firstneigh_half[i], + &numneigh_full[i],firstneigh_full[i], + &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], + rho0,&arho1[0][0],&arho2[0][0],arho2b, + &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0], + &errorflag); + if (errorflag) { + char str[128]; + sprintf(str,"MEAM library error %d",errorflag); + error->one(str); + } + offset += numneigh_half[i]; + } + + comm->reverse_comm_pair(this); + + meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, + &eng_vdwl,eatom,&ntype,type,fmap, + &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], + &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1, + dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); + if (errorflag) { + char str[128]; + sprintf(str,"MEAM library error %d",errorflag); + error->one(str); + } + + comm->comm_pair(this); + + offset = 0; + + // vptr is first value in vatom if it will be used by meam_force() + // else vatom may not exist, so pass dummy ptr + + double *vptr; + if (vflag_atom) vptr = &vatom[0][0]; + else vptr = &cutmax; + + for (ii = 0; ii < inum_half; ii++) { + i = ilist_half[ii]; + ifort = i+1; + meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, + &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], + &numneigh_half[i],firstneigh_half[i], + &numneigh_full[i],firstneigh_full[i], + &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], + dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, + &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], + &t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag); + if (errorflag) { + char str[128]; + sprintf(str,"MEAM library error %d",errorflag); + error->one(str); + } + offset += numneigh_half[i]; + } + + // change neighbor list indices back to C indexing + + neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); + neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); + + if (vflag_fdotr) virial_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairMEAM::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + map = new int[n+1]; + fmap = new int[n]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairMEAM::settings(int narg, char **arg) +{ + if (narg != 0) error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairMEAM::coeff(int narg, char **arg) +{ + int i,j,m,n; + + if (!allocated) allocate(); + + if (narg < 6) error->all("Incorrect args for pair coefficients"); + + // insure I,J args are * * + + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all("Incorrect args for pair coefficients"); + + // read MEAM element names between 2 filenames + // nelements = # of MEAM elements + // elements = list of unique element names + + if (nelements) { + for (i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + delete [] mass; + } + nelements = narg - 4 - atom->ntypes; + if (nelements < 1) error->all("Incorrect args for pair coefficients"); + elements = new char*[nelements]; + mass = new double[nelements]; + + for (i = 0; i < nelements; i++) { + n = strlen(arg[i+3]) + 1; + elements[i] = new char[n]; + strcpy(elements[i],arg[i+3]); + } + + // read MEAM library and parameter files + // pass all parameters to MEAM package + // tell MEAM package that setup is done + + read_files(arg[2],arg[2+nelements+1]); + meam_setup_done_(&cutmax); + + // read args that map atom types to MEAM elements + // map[i] = which element the Ith atom type is, -1 if not mapped + + for (i = 4 + nelements; i < narg; i++) { + m = i - (4+nelements) + 1; + for (j = 0; j < nelements; j++) + if (strcmp(arg[i],elements[j]) == 0) break; + if (j < nelements) map[m] = j; + else if (strcmp(arg[i],"NULL") == 0) map[m] = -1; + else error->all("Incorrect args for pair coefficients"); + } + + // clear setflag since coeff() called once with I,J = * * + + n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + // set mass for i,i in atom class + + int count = 0; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + if (i == j) atom->set_mass(i,mass[map[i]]); + count++; + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairMEAM::init_style() +{ + if (force->newton_pair == 0) + error->all("Pair style MEAM requires newton pair on"); + + // need full and half neighbor list + + int irequest_full = neighbor->request(this); + neighbor->requests[irequest_full]->id = 1; + neighbor->requests[irequest_full]->half = 0; + neighbor->requests[irequest_full]->full = 1; + int irequest_half = neighbor->request(this); + neighbor->requests[irequest_half]->id = 2; + neighbor->requests[irequest_half]->half = 0; + neighbor->requests[irequest_half]->half_from_full = 1; + neighbor->requests[irequest_half]->otherlist = irequest_full; + + // setup Fortran-style mapping array needed by MEAM package + // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index + // if type I is not a MEAM atom, fmap stores a 0 + + for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1; +} + +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + half or full +------------------------------------------------------------------------- */ + +void PairMEAM::init_list(int id, NeighList *ptr) +{ + if (id == 1) listfull = ptr; + else if (id == 2) listhalf = ptr; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairMEAM::init_one(int i, int j) +{ + return cutmax; +} + +/* ---------------------------------------------------------------------- */ + +void PairMEAM::read_files(char *globalfile, char *userfile) +{ + // open global meamf file on proc 0 + + FILE *fp; + if (comm->me == 0) { + fp = fopen(globalfile,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open MEAM potential file %s",globalfile); + error->one(str); + } + } + + // allocate parameter arrays + + int params_per_line = 19; + + int *lat = new int[nelements]; + int *ielement = new int[nelements]; + int *ibar = new int[nelements]; + double *z = new double[nelements]; + double *atwt = new double[nelements]; + double *alpha = new double[nelements]; + double *b0 = new double[nelements]; + double *b1 = new double[nelements]; + double *b2 = new double[nelements]; + double *b3 = new double[nelements]; + double *alat = new double[nelements]; + double *esub = new double[nelements]; + double *asub = new double[nelements]; + double *t0 = new double[nelements]; + double *t1 = new double[nelements]; + double *t2 = new double[nelements]; + double *t3 = new double[nelements]; + double *rozero = new double[nelements]; + + bool *found = new bool[nelements]; + for (int i = 0; i < nelements; i++) found[i] = false; + + // read each set of params from global MEAM file + // one set of params can span multiple lines + // store params if element name is in element list + // if element name appears multiple times, only store 1st entry + + int i,n,nwords; + char **words = new char*[params_per_line+1]; + char line[MAXLINE],*ptr; + int eof = 0; + + int nset = 0; + while (1) { + if (comm->me == 0) { + ptr = fgets(line,MAXLINE,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + + // strip comment, skip line if blank + + if (ptr = strchr(line,'#')) *ptr = '\0'; + nwords = atom->count_words(line); + if (nwords == 0) continue; + + // concatenate additional lines until have params_per_line words + + while (nwords < params_per_line) { + n = strlen(line); + if (comm->me == 0) { + ptr = fgets(&line[n],MAXLINE-n,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + if (ptr = strchr(line,'#')) *ptr = '\0'; + nwords = atom->count_words(line); + } + + if (nwords != params_per_line) + error->all("Incorrect format in MEAM potential file"); + + // words = ptrs to all words in line + // strip single and double quotes from words + + nwords = 0; + words[nwords++] = strtok(line,"' \t\n\r\f"); + while (words[nwords++] = strtok(NULL,"' \t\n\r\f")) continue; + + // skip if element name isn't in element list + + for (i = 0; i < nelements; i++) + if (strcmp(words[0],elements[i]) == 0) break; + if (i == nelements) continue; + + // skip if element already appeared + + if (found[i] == true) continue; + found[i] = true; + + // map lat string to an integer + + if (strcmp(words[1],"fcc") == 0) lat[i] = FCC; + else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC; + else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP; + else if (strcmp(words[1],"dim") == 0) lat[i] = DIM; + else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND; + else error->all("Unrecognized lattice type in MEAM file 1"); + + // store parameters + + z[i] = atof(words[2]); + ielement[i] = atoi(words[3]); + atwt[i] = atof(words[4]); + alpha[i] = atof(words[5]); + b0[i] = atof(words[6]); + b1[i] = atof(words[7]); + b2[i] = atof(words[8]); + b3[i] = atof(words[9]); + alat[i] = atof(words[10]); + esub[i] = atof(words[11]); + asub[i] = atof(words[12]); + t0[i] = atof(words[13]); + t1[i] = atof(words[14]); + t2[i] = atof(words[15]); + t3[i] = atof(words[16]); + rozero[i] = atof(words[17]); + ibar[i] = atoi(words[18]); + + nset++; + } + + // error if didn't find all elements in file + + if (nset != nelements) + error->all("Did not find all elements in MEAM library file"); + + // pass element parameters to MEAM package + + meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, + alat,esub,asub,t0,t1,t2,t3,rozero,ibar); + + // set element masses + + for (i = 0; i < nelements; i++) mass[i] = atwt[i]; + + // clean-up memory + + delete [] words; + + delete [] lat; + delete [] ielement; + delete [] ibar; + delete [] z; + delete [] atwt; + delete [] alpha; + delete [] b0; + delete [] b1; + delete [] b2; + delete [] b3; + delete [] alat; + delete [] esub; + delete [] asub; + delete [] t0; + delete [] t1; + delete [] t2; + delete [] t3; + delete [] rozero; + delete [] found; + + // done if user param file is NULL + + if (strcmp(userfile,"NULL") == 0) return; + + // open user param file on proc 0 + + if (comm->me == 0) { + fp = fopen(userfile,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open MEAM potential file %s",userfile); + error->one(str); + } + } + + // read settings + // pass them one at a time to MEAM package + // match strings to list of corresponding ints + + int which; + double value; + int nindex,index[3]; + int maxparams = 6; + char **params = new char*[maxparams]; + int nparams; + + eof = 0; + while (1) { + if (comm->me == 0) { + ptr = fgets(line,MAXLINE,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + + // strip comment, skip line if blank + + if (ptr = strchr(line,'#')) *ptr = '\0'; + nparams = atom->count_words(line); + if (nparams == 0) continue; + + // words = ptrs to all words in line + + nparams = 0; + params[nparams++] = strtok(line,"=(), '\t\n\r\f"); + while (nparams < maxparams && + (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f"))) + continue; + nparams--; + + for (which = 0; which < nkeywords; which++) + if (strcmp(params[0],keywords[which]) == 0) break; + if (which == nkeywords) { + char str[128]; + sprintf(str,"Keyword %s in MEAM parameter file not recognized", + params[0]); + error->all(str); + } + nindex = nparams - 2; + for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]); + + // map lattce_meam value to an integer + + if (which == 4) { + if (strcmp(params[nparams-1],"fcc") == 0) value = FCC; + else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC; + else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP; + else if (strcmp(params[nparams-1],"dim") == 0) value = DIM; + else if (strcmp(params[nparams-1],"dia") == 0) value = DIAMOND; + else if (strcmp(params[nparams-1],"b1") == 0) value = B1; + else if (strcmp(params[nparams-1],"c11") == 0) value = C11; + else if (strcmp(params[nparams-1],"l12") == 0) value = L12; + else error->all("Unrecognized lattice type in MEAM file 2"); + } + else value = atof(params[nparams-1]); + + // pass single setting to MEAM package + + int errorflag = 0; + meam_setup_param_(&which,&value,&nindex,index,&errorflag); + if (errorflag) { + char str[128]; + sprintf(str,"MEAM library error %d",errorflag); + error->all(str); + } + } + + delete [] params; +} + +/* ---------------------------------------------------------------------- */ + +int PairMEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) +{ + int i,j,k,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = rho0[j]; + buf[m++] = rho1[j]; + buf[m++] = rho2[j]; + buf[m++] = rho3[j]; + buf[m++] = frhop[j]; + buf[m++] = gamma[j]; + buf[m++] = dgamma1[j]; + buf[m++] = dgamma2[j]; + buf[m++] = dgamma3[j]; + buf[m++] = arho2b[j]; + buf[m++] = arho1[j][0]; + buf[m++] = arho1[j][1]; + buf[m++] = arho1[j][2]; + buf[m++] = arho2[j][0]; + buf[m++] = arho2[j][1]; + buf[m++] = arho2[j][2]; + buf[m++] = arho2[j][3]; + buf[m++] = arho2[j][4]; + buf[m++] = arho2[j][5]; + for (k = 0; k < 10; k++) buf[m++] = arho3[j][k]; + buf[m++] = arho3b[j][0]; + buf[m++] = arho3b[j][1]; + buf[m++] = arho3b[j][2]; + buf[m++] = t_ave[j][0]; + buf[m++] = t_ave[j][1]; + buf[m++] = t_ave[j][2]; + buf[m++] = tsq_ave[j][0]; + buf[m++] = tsq_ave[j][1]; + buf[m++] = tsq_ave[j][2]; + } + return 38; +} + +/* ---------------------------------------------------------------------- */ + +void PairMEAM::unpack_comm(int n, int first, double *buf) +{ + int i,k,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + rho0[i] = buf[m++]; + rho1[i] = buf[m++]; + rho2[i] = buf[m++]; + rho3[i] = buf[m++]; + frhop[i] = buf[m++]; + gamma[i] = buf[m++]; + dgamma1[i] = buf[m++]; + dgamma2[i] = buf[m++]; + dgamma3[i] = buf[m++]; + arho2b[i] = buf[m++]; + arho1[i][0] = buf[m++]; + arho1[i][1] = buf[m++]; + arho1[i][2] = buf[m++]; + arho2[i][0] = buf[m++]; + arho2[i][1] = buf[m++]; + arho2[i][2] = buf[m++]; + arho2[i][3] = buf[m++]; + arho2[i][4] = buf[m++]; + arho2[i][5] = buf[m++]; + for (k = 0; k < 10; k++) arho3[i][k] = buf[m++]; + arho3b[i][0] = buf[m++]; + arho3b[i][1] = buf[m++]; + arho3b[i][2] = buf[m++]; + t_ave[i][0] = buf[m++]; + t_ave[i][1] = buf[m++]; + t_ave[i][2] = buf[m++]; + tsq_ave[i][0] = buf[m++]; + tsq_ave[i][1] = buf[m++]; + tsq_ave[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int PairMEAM::pack_reverse_comm(int n, int first, double *buf) +{ + int i,k,m,last,size; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = rho0[i]; + buf[m++] = arho2b[i]; + buf[m++] = arho1[i][0]; + buf[m++] = arho1[i][1]; + buf[m++] = arho1[i][2]; + buf[m++] = arho2[i][0]; + buf[m++] = arho2[i][1]; + buf[m++] = arho2[i][2]; + buf[m++] = arho2[i][3]; + buf[m++] = arho2[i][4]; + buf[m++] = arho2[i][5]; + for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; + buf[m++] = arho3b[i][0]; + buf[m++] = arho3b[i][1]; + buf[m++] = arho3b[i][2]; + buf[m++] = t_ave[i][0]; + buf[m++] = t_ave[i][1]; + buf[m++] = t_ave[i][2]; + buf[m++] = tsq_ave[i][0]; + buf[m++] = tsq_ave[i][1]; + buf[m++] = tsq_ave[i][2]; + } + + return 30; +} + +/* ---------------------------------------------------------------------- */ + +void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,k,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + rho0[j] += buf[m++]; + arho2b[j] += buf[m++]; + arho1[j][0] += buf[m++]; + arho1[j][1] += buf[m++]; + arho1[j][2] += buf[m++]; + arho2[j][0] += buf[m++]; + arho2[j][1] += buf[m++]; + arho2[j][2] += buf[m++]; + arho2[j][3] += buf[m++]; + arho2[j][4] += buf[m++]; + arho2[j][5] += buf[m++]; + for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; + arho3b[j][0] += buf[m++]; + arho3b[j][1] += buf[m++]; + arho3b[j][2] += buf[m++]; + t_ave[j][0] += buf[m++]; + t_ave[j][1] += buf[m++]; + t_ave[j][2] += buf[m++]; + tsq_ave[j][0] += buf[m++]; + tsq_ave[j][1] += buf[m++]; + tsq_ave[j][2] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairMEAM::memory_usage() +{ + double bytes = 11 * nmax * sizeof(double); + bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double); + bytes += 3 * maxneigh * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + toggle neighbor list indices between zero- and one-based values + needed for access by MEAM Fortran library +------------------------------------------------------------------------- */ + +void PairMEAM::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh) +{ + int i,j,ii,jnum; + int *jlist; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + for (j = 0; j < jnum; j++) jlist[j]--; + } +} + +void PairMEAM::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh) +{ + int i,j,ii,jnum; + int *jlist; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + for (j = 0; j < jnum; j++) jlist[j]++; + } +} diff -Naur lammps-18Aug09/src/MEAM/pair_meam.h lammps-19Aug09/src/MEAM/pair_meam.h --- lammps-18Aug09/src/MEAM/pair_meam.h 2008-02-19 09:23:58.000000000 -0700 +++ lammps-19Aug09/src/MEAM/pair_meam.h 2009-08-18 11:01:49.000000000 -0600 @@ -1,96 +1,97 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifndef PAIR_MEAM_H -#define PAIR_MEAM_H - -extern "C" { - void meam_setup_global_(int *, int *, double *, int *, double *, double *, - double *, double *, double *, double *, double *, - double *, double *, double *, double *, double *, - double *, double *, int *); - void meam_setup_param_(int *, double *, int *, int *, int *); - void meam_setup_done_(double *); - - void meam_dens_init_(int *, int *, int *, int *, int *, - double *, int *, int *, int *, int *, - double *, double *, double *, double *, - double *, double *, - double *, double *, double *, double *, int *); - - void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, - int *, int *, int *, - double *, double *, double *, double *, - double *, double *, - double *, double *, double *, double *, - double *, double *, - double *, double *, double *, double *, int *); - - void meam_force_(int *, int *, int *, int *, int *, int *, - double *, double *, int *, int *, int *, - double *, int *, int *, int *, int *, double *, double *, - double *, double *, double *, double *, double *, double *, - double *, double *, double *, double *, double *, double *, - double *, double *, double *, double *, double *, int *); - - void meam_cleanup_(); -} - - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairMEAM : public Pair { - public: - PairMEAM(class LAMMPS *); - ~PairMEAM(); - void compute(int, int); - void settings(int, char **); - void coeff(int, char **); - void init_style(); - void init_list(int, class NeighList *); - double init_one(int, int); - - int pack_comm(int, int *, double *, int, int *); - void unpack_comm(int, int, double *); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); - double memory_usage(); - - private: - double cutmax; // max cutoff for all elements - int nelements; // # of unique elements - char **elements; // names of unique elements - double *mass; // mass of each element - - int *map; // mapping from atom types to elements - int *fmap; // Fortran version of map array for MEAM lib - - int maxneigh; - double *scrfcn,*dscrfcn,*fcpair; - - int nmax; - double *rho,*rho0,*rho1,*rho2,*rho3,*frhop; - double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b; - double **arho1,**arho2,**arho3,**arho3b,**t_ave; - - void allocate(); - void read_files(char *, char *); - void neigh_f2c(int, int *, int *, int **); - void neigh_c2f(int, int *, int *, int **); -}; - -} - -#endif +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef PAIR_MEAM_H +#define PAIR_MEAM_H + +extern "C" { + void meam_setup_global_(int *, int *, double *, int *, double *, double *, + double *, double *, double *, double *, double *, + double *, double *, double *, double *, double *, + double *, double *, int *); + void meam_setup_param_(int *, double *, int *, int *, int *); + void meam_setup_done_(double *); + + void meam_dens_init_(int *, int *, int *, int *, int *, + double *, int *, int *, int *, int *, + double *, double *, double *, double *, + double *, double *, + double *, double *, double *, double *, double *, + int *); + + void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, + int *, int *, int *, + double *, double *, double *, double *, + double *, double *, double *, + double *, double *, double *, double *, + double *, double *, + double *, double *, double *, double *, int *); + + void meam_force_(int *, int *, int *, int *, int *, int *, + double *, double *, int *, int *, int *, + double *, int *, int *, int *, int *, double *, double *, + double *, double *, double *, double *, double *, double *, + double *, double *, double *, double *, double *, double *, + double *, double *, double *, double *, double *, double *, int *); + + void meam_cleanup_(); +} + + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairMEAM : public Pair { + public: + PairMEAM(class LAMMPS *); + ~PairMEAM(); + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + void init_style(); + void init_list(int, class NeighList *); + double init_one(int, int); + + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + double memory_usage(); + + private: + double cutmax; // max cutoff for all elements + int nelements; // # of unique elements + char **elements; // names of unique elements + double *mass; // mass of each element + + int *map; // mapping from atom types to elements + int *fmap; // Fortran version of map array for MEAM lib + + int maxneigh; + double *scrfcn,*dscrfcn,*fcpair; + + int nmax; + double *rho,*rho0,*rho1,*rho2,*rho3,*frhop; + double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b; + double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave; + + void allocate(); + void read_files(char *, char *); + void neigh_f2c(int, int *, int *, int **); + void neigh_c2f(int, int *, int *, int **); +}; + +} + +#endif diff -Naur lammps-18Aug09/src/min_linesearch.cpp lammps-19Aug09/src/min_linesearch.cpp --- lammps-18Aug09/src/min_linesearch.cpp 2009-08-18 10:13:16.000000000 -0600 +++ lammps-19Aug09/src/min_linesearch.cpp 2009-08-18 11:13:07.000000000 -0600 @@ -248,26 +248,7 @@ // backtrack with alpha until energy decrease is sufficient while (1) { - if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - x0atom = x0extra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] = x0atom[i]; - } - if (nextra_global) modify->min_step(alpha,hextra); - for (i = 0; i < n3; i++) x[i] += alpha*h[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - hatom = hextra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i]; - } - ecurrent = energy_force(1); - nfunc++; + ecurrent = alpha_step(alpha,1,nfunc); // if energy change is better than ideal, exit with success @@ -283,17 +264,7 @@ // reset to starting point, exit with error if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) { - if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - x0atom = x0extra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] = x0atom[i]; - } - ecurrent = energy_force(0); - nfunc++; + ecurrent = alpha_step(0.0,0,nfunc); return ZEROALPHA; } } @@ -398,26 +369,7 @@ alphaprev = 0.0; while (1) { - if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - x0atom = x0extra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] = x0atom[i]; - } - if (nextra_global) modify->min_step(alpha,hextra); - for (i = 0; i < n3; i++) x[i] += alpha*h[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - hatom = hextra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i]; - } - ecurrent = energy_force(1); - nfunc++; + ecurrent = alpha_step(alpha,1,nfunc); // compute new fh, alpha, delfh @@ -455,17 +407,7 @@ // if fh or delfh is epsilon, reset to starting point, exit with error if (fabs(fh) < EPS_QUAD || fabs(delfh) < EPS_QUAD) { - if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - x0atom = x0extra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] = x0atom[i]; - } - ecurrent = energy_force(0); - nfunc++; + ecurrent = alpha_step(0.0,0,nfunc); return ZEROQUAD; } @@ -477,27 +419,7 @@ alpha0 = alpha - (alpha-alphaprev)*fh/delfh; if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) { - if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - x0atom = x0extra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] = x0atom[i]; - } - if (nextra_global) modify->min_step(alpha0,hextra); - for (i = 0; i < n3; i++) x[i] += alpha0*h[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - hatom = hextra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] += alpha0*hatom[i]; - } - ecurrent = energy_force(1); - nfunc++; - + ecurrent = alpha_step(alpha0,1,nfunc); return 0; } @@ -521,18 +443,47 @@ // reset to starting point, exit with error if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) { - if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; - if (nextra_atom) - for (m = 0; m < nextra_atom; m++) { - xatom = xextra_atom[m]; - x0atom = x0extra_atom[m]; - n = extra_nlen[m]; - for (i = 0; i < n; i++) xatom[i] = x0atom[i]; - } - ecurrent = energy_force(0); - nfunc++; + ecurrent = alpha_step(0.0,0,nfunc); return ZEROALPHA; } } } + +/* ---------------------------------------------------------------------- */ + +double MinLineSearch::alpha_step(double alpha, int resetflag, int &nfunc) +{ + int i,n,m; + double *xatom,*x0atom,*hatom; + + // reset to starting point + + if (nextra_global) modify->min_step(0.0,hextra); + for (i = 0; i < n3; i++) x[i] = x0[i]; + if (nextra_atom) + for (m = 0; m < nextra_atom; m++) { + xatom = xextra_atom[m]; + x0atom = x0extra_atom[m]; + n = extra_nlen[m]; + for (i = 0; i < n; i++) xatom[i] = x0atom[i]; + } + + // step forward along h + + if (alpha > 0.0) { + if (nextra_global) modify->min_step(alpha,hextra); + for (i = 0; i < n3; i++) x[i] += alpha*h[i]; + if (nextra_atom) + for (m = 0; m < nextra_atom; m++) { + xatom = xextra_atom[m]; + hatom = hextra_atom[m]; + n = extra_nlen[m]; + for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i]; + } + } + + // compute and return new energy + + nfunc++; + return energy_force(resetflag); +} diff -Naur lammps-18Aug09/src/min_linesearch.h lammps-19Aug09/src/min_linesearch.h --- lammps-18Aug09/src/min_linesearch.h 2009-08-14 14:54:10.000000000 -0600 +++ lammps-19Aug09/src/min_linesearch.h 2009-08-18 11:13:07.000000000 -0600 @@ -46,6 +46,8 @@ FnPtr linemin; int linemin_backtrack(double, double &, int &); int linemin_quadratic(double, double &, int &); + + double alpha_step(double, int, int &); }; }