添え字いぱーい - ita’s diary

ちょいと計算しますよ。force of MEAM potential
$E= F_s(\rho_s)+F_p(\rho_p)+F_d(\rho_d)$
$\rho_s = \sum_j \phi_s(r_j)$
$\rho_p = \sqrt{\rho_{\nu}\rho_{\nu}},\,\,\, \rho_\nu = \sum_j \phi_p(r_j) \nu_j/r_j$
$\rho_d = \sqrt{\rho_{\nu\mu}\rho_{\nu\mu}},\,\,\, \rho_{\nu\mu} = \sum_j \phi_d(r_j) (\nu_j\mu_j/r_j^2 -\delta(\nu,\mu)/3)$
$r_j = \sqrt{\nu_j \nu_j}$
$\partial_\lambda F_s = F^\prime_s(\rho_s) \partial_\lambda \rho_s, \,\, \partial_\lambda \rho_s=\sum_j \partial_\lambda \phi_s(r_j) = \sum_j \phi^\prime_s(r_j) \lambda_j/r_j$
$\partial_\lambda F_p = F^\prime_p(\rho_p) \partial_\lambda \rho_p,\,\,\, \partial_\lambda \rho_p =(\rho_\nu \rho_{\nu,\lambda})/\rho_p$
$\rho_{\nu,\lambda}=\sum_j (\phi_p^\prime(r_j) -\phi_p(r_j)/r_j)\frac{\nu_j\lambda_j}{r^2_j} + \delta(\nu,\lambda)\phi_p(r_j)/r_j$

$\partial_\lambda \rho_d =(\rho_{\nu\mu} \rho_{\nu\mu,\lambda})/\rho_d$
$\rho_{\nu\mu,\lambda}=\sum_j \phi^\prime_d(r_j) \lambda_j/r_j(\nu_j\mu_j/r_j^2 -\delta(\nu,mu)/3) +\phi_d(r_j) (\delta(\lambda,nu)\mu_j r_j^2+\delta(\lambda,\mu)\nu_j r^2_j -2\nu_j\mu_j \lambda_j)/r_j^4$

// electron density function
extern double PhiS(double r);
extern double PhiP(double r);
extern double PhiD(double r);
// and their delivative w.r.t. r
extern double DPhiS(double r);
extern double DPhiP(double r);
extern double DPhiD(double r);

class Atom
{
  public:
  double Pos[3];
  double Force[3];
  //electron density
  double RhoS, RhoP, RhoP_n[3], RhoD, RhoD_nm[3][3];
  // gradient
  double RhoS_l[3], RhoP_l[3], RhoP_nl[3][3], RhoD_l[3], RhoD_nml[3][3][3];

  double Fs, Fp, Fd, DFs, DFp, DFd;

};
class Pair
{
  public:
  class Atom *Atom1, *Atom2;
  double R1, R[3];

  void CalcRho()
  {
     int n,m,l;
     R1=0.0;
     for(n=0;n<3;n++)
     {
       R[n]=Atom1->Pos[n] - Atom2->Pos[n];
       R1 += R[n]*R[n];
     }
     R1 = sqrt(R1);
     for(n=0;n<3;n++) R[n] /= R1;

     // S-part
     double RhoS = PhiS(R1);
     double DRhoS = DPhiS(R1);
     Atom1->RhoS += RhoS;
     Atom2->RhoS += RhoS;
     for(l=0;l<3;l++)
     {
            Atom1->RhoS_l[l] += DRhoS*R[l];
            Atom2->RhoS_l[l] -= DRhoS*R[l];
     }

     // P-part
     double RhoP = PhiP(R1);
     double DRhoP = DPhiP(R1);
     for(n=0;n<3;n++)
     {
       Atom1->RhoP_n[n] += RhoP*R[n];
       Atom2->RhoP_n[n] -= RhoP*R[n];
       for(l=0;l<3;l++)
       {
          double tmp=(DRhoP-RhoP/R1)*R[n]*R[l]
          if (l==n) tmp += RhoP/R1;
          Atom1->RhoP_nl[n][l] += tmp;
          Atom2->RhoP_nl[n][l] -= tmp;
       }
     }//n

     double RhoD = PhiD(R1);
     double DRhoD = DPhiD(R1);
     for(n=0;n<3;n++) {
     for(m=0;m<3;m++) {
        double tmp1 = R[n] * R[m];
        if (n==m) tmp1 -= 1.0/3.0;
        Atom1->RhoD_nm[n][m] += tmp1*RhoD;
        Atom2->RhoD_nm[n][m] += tmp1*RhoD;

        for(l=0;l<3;l++)
        {
            double tmp2 = DRhoD * R[l] * tmp1;
            double tmp3 = -2*R[l]*R[n]*R[m];
            if (l==n) tmp3 += R[n];
            if (l==m) tmp3 += R[m];
            tmp2 += tmp3 * RhoD/R1;
            Atom1->RhoD_nml[n][m][l] += tmp2;
            Atom2->RhoD_nml[n][m][l] -= tmp2;
       } //l
     }} // n m
   }
};