use cube() and square() from math_special.h instead of pow(x,3.0) and pow(x,2.0) (71ce1c20) · Commits · 郑智淋 / lammps

src/MANYBODY/pair_tersoff.cpp

+5 −3

Original line number	Diff line number	Diff line
		@@ -30,9 +30,11 @@
		#include "error.h"

		#include "math_const.h"
		#include "math_special.h"

		using namespace LAMMPS_NS;
		using namespace MathConst;
		using namespace MathSpecial;

		#define MAXLINE 1024
		#define DELTA 4
		@@ -604,7 +606,7 @@ double PairTersoff::zeta(Param *param, double rsqij, double rsqik,
		costheta = (delrij[0]delrik[0] + delrij[1]delrik[1] +
		delrij[2]delrik[2]) / (rijrik);

		if (param->powermint == 3) arg = pow(param->lam3 * (rij-rik),3.0);
		if (param->powermint == 3) arg = cube(param->lam3 * (rij-rik));
		else arg = param->lam3 * (rij-rik);

		if (arg > 69.0776) ex_delr = 1.e30;
		@@ -744,7 +746,7 @@ void PairTersoff::ters_zetaterm_d(double prefactor,

		fc = ters_fc(rik,param);
		dfc = ters_fc_d(rik,param);
		if (param->powermint == 3) tmp = pow(param->lam3 * (rij-rik),3.0);
		if (param->powermint == 3) tmp = cube(param->lam3 * (rij-rik));
		else tmp = param->lam3 * (rij-rik);

		if (tmp > 69.0776) ex_delr = 1.e30;
		@@ -752,7 +754,7 @@ void PairTersoff::ters_zetaterm_d(double prefactor,
		else ex_delr = exp(tmp);

		if (param->powermint == 3)
		ex_delr_d = 3.0pow(param->lam3,3.0) pow(rij-rik,2.0)*ex_delr;
		ex_delr_d = 3.0cube(param->lam3) square(rij-rik)*ex_delr;
		else ex_delr_d = param->lam3 * ex_delr;

		cos_theta = vec3_dot(rij_hat,rik_hat);

+5 −3

Original line number	Diff line number	Diff line
		@@ -25,11 +25,13 @@
		#include "force.h"
		#include "comm.h"
		#include "math_const.h"
		#include "math_special.h"
		#include "memory.h"
		#include "error.h"

		using namespace LAMMPS_NS;
		using namespace MathConst;
		using namespace MathSpecial;

		#define MAXLINE 1024
		#define DELTA 4
		@@ -241,7 +243,7 @@ double PairTersoffMOD::zeta(Param *param, double rsqij, double rsqik,
		costheta = (delrij[0]delrik[0] + delrij[1]delrik[1] +
		delrij[2]delrik[2]) / (rijrik);

		if (param->powermint == 3) arg = pow(param->lam3 * (rij-rik),3.0);
		if (param->powermint == 3) arg = cube(param->lam3 * (rij-rik));
		else arg = param->lam3 * (rij-rik);

		if (arg > 69.0776) ex_delr = 1.e30;
		@@ -314,7 +316,7 @@ void PairTersoffMOD::ters_zetaterm_d(double prefactor,

		fc = ters_fc(rik,param);
		dfc = ters_fc_d(rik,param);
		if (param->powermint == 3) tmp = pow(param->lam3 * (rij-rik),3.0);
		if (param->powermint == 3) tmp = cube(param->lam3 * (rij-rik));
		else tmp = param->lam3 * (rij-rik);

		if (tmp > 69.0776) ex_delr = 1.e30;
		@@ -322,7 +324,7 @@ void PairTersoffMOD::ters_zetaterm_d(double prefactor,
		else ex_delr = exp(tmp);

		if (param->powermint == 3)
		ex_delr_d = 3.0pow(param->lam3,3.0) pow(rij-rik,2.0)*ex_delr;
		ex_delr_d = 3.0cube(param->lam3) square(rij-rik)*ex_delr;
		else ex_delr_d = param->lam3 * ex_delr;

		cos_theta = vec3_dot(rij_hat,rik_hat);

+4 −2

Original line number	Diff line number	Diff line
		@@ -28,9 +28,11 @@
		#include "memory.h"
		#include "error.h"
		#include "math_const.h"
		#include "math_special.h"

		using namespace LAMMPS_NS;
		using namespace MathConst;
		using namespace MathSpecial;

		#define MAXLINE 1024
		#define DELTA 4
		@@ -226,7 +228,7 @@ void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce,

		// ZBL repulsive portion

		double esq = pow(global_e,2.0);
		double esq = square(global_e);
		double a_ij = (0.8854*global_a_0) /
		(pow(param->Z_i,0.23) + pow(param->Z_j,0.23));
		double premult = (param->Z_i * param->Z_j * esq)/(4.0MY_PIglobal_epsilon_0);
		@@ -286,5 +288,5 @@ double PairTersoffZBL::F_fermi(double r, Param *param)
		double PairTersoffZBL::F_fermi_d(double r, Param *param)
		{
		return param->ZBLexpscaleexp(-param->ZBLexpscale(r-param->ZBLcut)) /
		pow(1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)),2.0);
		square(1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)));
		}