--- forces-edip.1.2.f Fri Apr 19 10:42:33 2002 +++ forces-edip.1.3.f Fri Apr 19 10:42:49 2002 @@ -2,7 +2,7 @@ C forces-edip.f C ------------- -C Version 1.2f +C Version 1.3f C Force and Energy Calculation with the C Environment-Dependent Interatomic Potential @@ -123,7 +123,7 @@ par_bg=par_a par_eta = par_delta/par_Qo pot_cutoff = par_a - delta_safe = 0.2 + delta_safe = 0.2d0 end subroutine @@ -136,21 +136,21 @@ print *, "EDIP-Si parameters:" C taken from Justo et al., Phys. Rev. B 58, 2539 (1998). - par_cap_A = 5.6714030 - par_cap_B = 2.0002804 - par_rh = 1.2085196 - par_a = 3.1213820 - par_sig = 0.5774108 - par_lam = 1.4533108 - par_gam = 1.1247945 - par_b = 3.1213820 - par_c = 2.5609104 - par_delta = 78.7590539 - par_mu = -0.6966326 - par_Qo = 312.1341346 - par_palp = 1.4074424 - par_bet = -0.0070975 - par_alp = 3.1083847 + par_cap_A = 5.6714030d0 + par_cap_B = 2.0002804d0 + par_rh = 1.2085196d0 + par_a = 3.1213820d0 + par_sig = 0.5774108d0 + par_lam = 1.4533108d0 + par_gam = 1.1247945d0 + par_b = 3.1213820d0 + par_c = 2.5609104d0 + par_delta = 78.7590539d0 + par_mu = -0.6966326d0 + par_Qo = 312.1341346d0 + par_palp = 1.4074424d0 + par_bet = -0.0070975d0 + par_alp = 3.1083847d0 print *, par_cap_A,par_cap_B,par_rh,par_a,par_sig print *, par_lam,par_gam,par_b,par_c,par_delta print *, par_mu,par_Qo,par_palp,par_bet,par_alp @@ -244,18 +244,18 @@ L_z_div_2 = L_z/2.0D0 do i=1, N_own - f(1,i) = 0.0 - f(2,i) = 0.0 - f(3,i) = 0.0 + f(1,i) = 0.0d0 + f(2,i) = 0.0d0 + f(3,i) = 0.0d0 end do - coord_total=0.0 - virial=0.0 - virial_xyz(1)= 0.0 - virial_xyz(2)= 0.0 - virial_xyz(3)= 0.0 - V2=0.0 - V3=0.0 + coord_total=0.0d0 + virial=0.0d0 + virial_xyz(1)= 0.0d0 + virial_xyz(2)= 0.0d0 + virial_xyz(3)= 0.0d0 + V2=0.0d0 + V3=0.0d0 C COMBINE COEFFICIENTS @@ -275,7 +275,7 @@ C RESET COORDINATION AND NEIGHBOR NUMBERS - Z = 0.0 + Z = 0.0d0 n2 = 0 n3 = 0 nz = 0 @@ -316,7 +316,7 @@ rsqr = dx*dx + dy*dy + dz*dz if(rsqr < asqr) then r = sqrt(rsqr) - rinv = 1.0/r + rinv = 1.0d0/r dx = dx * rinv dy = dy * rinv dz = dz * rinv @@ -325,7 +325,7 @@ n2 = n2 + 1 num2(n2) = j - rmainv = 1.0/(r-par_a) + rmainv = 1.0d0/(r-par_a) s2_t0(n2) = par_cap_A*dexp(par_sig*rmainv) s2_t1(n2) = (par_cap_B*rinv)**par_rh s2_t2(n2) = par_rh*rinv @@ -342,7 +342,7 @@ n3 = n3 + 1 num3(n3) = j - rmbinv = 1.0/(r-par_bg) + rmbinv = 1.0d0/(r-par_bg) temp1 = par_gam*rmbinv temp0 = dexp(temp1) s3_g(n3) = temp0 @@ -360,13 +360,13 @@ if(r .lt. par_b) then if(r .lt. par_c) then - Z = Z + 1.0 + Z = Z + 1.0d0 else nz = nz + 1 numz(nz) = j xinv = bmc/(r-par_c) xinv3 = xinv*xinv*xinv - den = 1.0/(1 - xinv3) + den = 1.0d0/(1.d0 - xinv3) temp1 = par_alp*den fZ = dexp(temp1) Z = Z + fZ @@ -399,7 +399,7 @@ Z = tricks_Zfix coord_total = coord_total + Z pZ = par_palp*dexp(par_bet*Z*Z) - dp = 0.0 + dp = 0.0d0 else @@ -475,10 +475,10 @@ if(fixZ .ne. 0) then winv = Qort*dexp(muhalf*Z) - dwinv = 0.0 + dwinv = 0.0d0 temp0 = dexp(-u4*Z) tau = u1+u2*temp0*(u3-temp0) - dtau = 0.0 + dtau = 0.0d0 else C COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION @@ -490,7 +490,7 @@ temp0 = exp(-u4*Z) tau = u1+u2*temp0*(u3-temp0) C -cosine of angular minimum - dtau = u5*temp0*(2*temp0-u3) + dtau = u5*temp0*(2.d0*temp0-u3) C its derivative end if @@ -516,8 +516,8 @@ x = (lcos + tau)*winv temp0 = exp(-x*x) - H = par_lam*(1 - temp0 + par_eta*x*x) - dHdx = 2*par_lam*x*(temp0 + par_eta) + H = par_lam*(1.d0 - temp0 + par_eta*x*x) + dHdx = 2.d0*par_lam*x*(temp0 + par_eta) dhdl = dHdx*winv @@ -642,18 +642,18 @@ C dE/dZ*dZ/dr contribution to virial virial = virial - sz_r(nl) * (dEdrlx*sz_dx(nl)+ - & dEdrly*sz_dy(nl) + dEdrlz*sz_dz(nl)) + & dEdrly*sz_dy(nl) + dEdrlz*sz_dz(nl)) virial_xyz(1) = virial_xyz(1) - - & sz_r(nl)*(dEdrlx*sz_dx(nl)) + & sz_r(nl)*(dEdrlx*sz_dx(nl)) virial_xyz(2) = virial_xyz(2) - - & sz_r(nl)*(dEdrly*sz_dy(nl)) + & sz_r(nl)*(dEdrly*sz_dy(nl)) virial_xyz(3) = virial_xyz(3) - - & sz_r(nl)*(dEdrlz*sz_dz(nl)) + & sz_r(nl)*(dEdrlz*sz_dz(nl)) end do end if end do E_potential = V2+V3 - virial = virial / 3.0 + virial = virial / 3.0d0 end subroutine