--- forces-edip.1.11.f Fri Apr 19 10:44:38 2002 +++ forces-edip.1.12.f Fri Apr 19 10:45:34 2002 @@ -2,7 +2,7 @@ C forces-edip.f C ------------- -C Version 1.11f +C Version 1.12f C Force and Energy Calculation with the C Environment-Dependent Interatomic Potential @@ -110,22 +110,34 @@ C ****************************************** - subroutine compute_forces_EDIP(N_own, pos, L_x, L_y, L_z, - & E_potential, f) - implicit none + subroutine edip_casc + & (natms,engsrp,engtbp,virtbp,cell,list,lentry, + & xxx,yyy,zzz,fxx,fyy,fzz) C NOTE: you do not need these header files, since you will have to C customize the way that atoms are passed to the force subroutine C for your particular MD program. - include "edip_neighbors_include.h" + include "dl_params.inc" +c implicit none +c integer msatms,mxlist,mxatms +c parameter (msatms=25000) +c parameter (mxlist=72) +c parameter (mxatms=1728) c#define FIXZ #undef FIXZ c switch on/off option using fixed Z +#define REP_NUCL +c#undef REP_NUCL +c switch on/off option calling external nuclear replusion term + C ------------------------- EDIP PARAMETERS --------------------------- + double precision par_unit + parameter (par_unit = 9648.531d0) + double precision par_cap_A,par_cap_B,par_rh,par_a,par_sig double precision par_lam,par_gam,par_b,par_c,par_delta double precision par_mu,par_Qo,par_palp,par_bet,par_alp @@ -134,12 +146,12 @@ c double precision pot_cutoff,delta_safe C taken from Justo et al., Phys. Rev. B 58, 2539 (1998). - parameter (par_cap_A = 5.6714030d0) + parameter (par_cap_A = 5.6714030d0*par_unit) parameter (par_cap_B = 2.0002804d0) parameter (par_rh = 1.2085196d0) parameter (par_a = 3.1213820d0) parameter (par_sig = 0.5774108d0) - parameter (par_lam = 1.4533108d0) + parameter (par_lam = 1.4533108d0*par_unit) parameter (par_gam = 1.1247945d0) parameter (par_b = 3.1213820d0) parameter (par_c = 2.5609104d0) @@ -172,11 +184,13 @@ C ------------------------- VARIABLE DECLARATIONS ------------------------- - integer N_own,MAX_NBRS_1 - double precision pos(3,N_own) - double precision E_potential - double precision f(3,N_own) - double precision L_x, L_y, L_z + integer natms,MAX_NBRS_1 + integer lentry(msatms),list(msatms,mxlist) + double precision engsrp,engtbp,virtbp + double precision cell(9), L_x, L_y, L_z + double precision xxx(mxatms),yyy(mxatms),zzz(mxatms) + double precision fxx(mxatms),fyy(mxatms),fzz(mxatms) + parameter (MAX_NBRS_1 = mxlist) integer i,j,k,l,n double precision dx,dy,dz,r,rsqr @@ -232,6 +246,10 @@ parameter (tricks_Zfix = 5) #endif +#ifdef REP_NUCL + double precision Vnucl,Fnucl +#endif + integer nj,nk,nl C indices for the store arrays @@ -241,14 +259,17 @@ double precision coord_total + L_x = cell(1) + L_y = cell(5) + L_z = cell(9) L_x_div_2 = L_x/2.0D0 L_y_div_2 = L_y/2.0D0 L_z_div_2 = L_z/2.0D0 - do i=1, N_own - f(1,i) = 0.0d0 - f(2,i) = 0.0d0 - f(3,i) = 0.0d0 + do i=1, natms + fxx(i) = 0.0d0 + fyy(i) = 0.0d0 + fzz(i) = 0.0d0 end do coord_total=0.0d0 @@ -262,7 +283,7 @@ C --- LEVEL 1: OUTER LOOP OVER ATOMS --- - do i=1, N_own + do i=1, natms C RESET COORDINATION AND NEIGHBOR NUMBERS @@ -274,13 +295,13 @@ C --- LEVEL 2: LOOP PREPASS OVER PAIRS --- - do n=p_nbrs(i), p_nbrs(i+1)-1 - j = neighbors(n) + do n = 1, lentry(i) + j = list(i,n) C TEST IF WITHIN OUTER CUTOFF - dx = pos(1,j)-pos(1,i) + dx = xxx(j)-xxx(i) if (dx .gt. L_x_div_2) then dx = dx - L_x else if (dx .lt. -L_x_div_2) then @@ -288,7 +309,7 @@ end if C dx periodic if(dabs(dx) < par_a) then - dy = pos(2,j)-pos(2,i) + dy = yyy(j)-yyy(i) if (dy .gt. L_y_div_2) then dy = dy - L_y else if (dy .lt. -L_y_div_2) then @@ -296,7 +317,7 @@ end if C dy periodic if(dabs(dy) < par_a) then - dz = pos(3,j)-pos(3,i) + dz = zzz(j)-zzz(i) if (dz .gt. L_z_div_2) then dz = dz - L_z else if (dz .lt. -L_z_div_2) then @@ -419,17 +440,27 @@ C two-body forces dV2j = - s2_t0(nj)*s2_t1(nj)*s2_t2(nj) - temp0*s2_t3(nj) + +#ifdef REP_NUCL + + call rep_nucl(s2_r(nj),Vnucl,Fnucl) + + V2 = V2 + Vnucl + dV2j = dV2j + Fnucl + +#endif + C dV2/dr dV2ijx = dV2j * s2_dx(nj) dV2ijy = dV2j * s2_dy(nj) dV2ijz = dV2j * s2_dz(nj) - f(1,i) = f(1,i) + dV2ijx - f(2,i) = f(2,i) + dV2ijy - f(3,i) = f(3,i) + dV2ijz + fxx(i) = fxx(i) + dV2ijx + fyy(i) = fyy(i) + dV2ijy + fzz(i) = fzz(i) + dV2ijz j = num2(nj) - f(1,j) = f(1,j) - dV2ijx - f(2,j) = f(2,j) - dV2ijy - f(3,j) = f(3,j) - dV2ijz + fxx(j) = fxx(j) - dV2ijx + fyy(j) = fyy(j) - dV2ijy + fzz(j) = fzz(j) - dV2ijz C dV2/dr contribution to virial @@ -544,15 +575,15 @@ C apply radial + angular forces to i, j, k - f(1,j) = f(1,j) - fjx - f(2,j) = f(2,j) - fjy - f(3,j) = f(3,j) - fjz - f(1,k) = f(1,k) - fkx - f(2,k) = f(2,k) - fky - f(3,k) = f(3,k) - fkz - f(1,i) = f(1,i) + fjx + fkx - f(2,i) = f(2,i) + fjy + fky - f(3,i) = f(3,i) + fjz + fkz + fxx(j) = fxx(j) - fjx + fyy(j) = fyy(j) - fjy + fzz(j) = fzz(j) - fjz + fxx(k) = fxx(k) - fkx + fyy(k) = fyy(k) - fky + fzz(k) = fzz(k) - fkz + fxx(i) = fxx(i) + fjx + fkx + fyy(i) = fyy(i) + fjy + fky + fzz(i) = fzz(i) + fjz + fkz C dV3/dR contributions to virial @@ -586,13 +617,13 @@ dEdrlx = dEdrl * sz_dx(nl) dEdrly = dEdrl * sz_dy(nl) dEdrlz = dEdrl * sz_dz(nl) - f(1,i) = f(1,i) + dEdrlx - f(2,i) = f(2,i) + dEdrly - f(3,i) = f(3,i) + dEdrlz + fxx(i) = fxx(i) + dEdrlx + fyy(i) = fyy(i) + dEdrly + fzz(i) = fzz(i) + dEdrlz l = numz(nl) - f(1,l) = f(1,l) - dEdrlx - f(2,l) = f(2,l) - dEdrly - f(3,l) = f(3,l) - dEdrlz + fxx(l) = fxx(l) - dEdrlx + fyy(l) = fyy(l) - dEdrly + fzz(l) = fzz(l) - dEdrlz C dE/dZ*dZ/dr contribution to virial @@ -610,6 +641,11 @@ end do - E_potential = V2+V3 - virial = virial / 3.0d0 +c E_potential = V2+V3 +c virial = virial / 3.0d0 + + engsrp = V2 + engtbp = V3 + virtbp = - virial + end subroutine