diff --git a/src/pair_zbl.cpp b/src/pair_zbl.cpp index a9882f1fe7..c44c9542b4 100644 --- a/src/pair_zbl.cpp +++ b/src/pair_zbl.cpp @@ -240,6 +240,10 @@ void PairZBL::init_style() cut_innersq = cut_inner * cut_inner; cut_globalsq = cut_global * cut_global; + // metal units: e^2/angstrom = 2*Ryd*bohr + // econv = 2.0*13.6058*0.529177 = 14.3998; + // all units: + econv = force->qqr2e; } /* ---------------------------------------------------------------------- @@ -273,9 +277,17 @@ double PairZBL::init_one(int i, int j) // sw1 = A // sw2 = B + // de2dr2 = 2*A*t + 3*B*t^2 + + // Require that at t = tc: + // e = -Fc + // dedr = -Fc' + // d2edr2 = -Fc'' + + // Hence: // A = (-3Fc' + tc*Fc'')/tc^2 // B = ( 2Fc' - tc*Fc'')/tc^3 - // C = = -Fc + tc/2*Fc' - tc^2/12*Fc'' + // C = -Fc + tc/2*Fc' - tc^2/12*Fc'' double tc = cut_global - cut_inner; double fc = e_zbl(cut_global, i, j); diff --git a/src/pair_zbl.h b/src/pair_zbl.h index 690ae8107d..c93106369f 100644 --- a/src/pair_zbl.h +++ b/src/pair_zbl.h @@ -37,7 +37,7 @@ class PairZBL : public Pair { protected: double cut_global, cut_inner; - double cut_globalsq, cut_innersq; + double cut_globalsq, cut_innersq, econv; double *z; double **d1a,**d2a,**d3a,**d4a,**zze; double **sw1,**sw2,**sw3,**sw4,**sw5; @@ -62,10 +62,6 @@ namespace PairZBLConstants { static const double d2 = 0.40290; static const double d3 = 0.94229; static const double d4 = 3.19980; - - // units: e^2/angstrom = 2*Ryd*bohr - - static const double econv = 2.0*13.6058*0.529177; } }