/* Fill in motion, pote, f */
void potential_force()
{
    register int i, j, k;
    int l, swapped;
    register double cc, ff;
    double h, sc, ss, tt, cu, cv;
    double ds[DIMENSION], dx[4], m[DIMENSION], xx[DIMENSION], wk[DIMENSION],
        wx[DIMENSION], force[DIMENSION], gorce[DIMENSION];
    motion = Config_analyze_motion (Config_Aapp_to_Alib, NULL);
    IZERO (np, nn);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            if (tp[i] == tp[j]) continue;
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( dx[3] < SQUARE(RCUT_ZrC) )
            {
                dx[3] = sqrt(dx[3]);
                if ((nn[i]==MAX_FIRST) || (nn[j]==MAX_FIRST))
                    pe("potential_force: MAX_FIRST=%d reached.\n", MAX_FIRST);
                V3DIV ( dx, dx[3], &NEIGH(i,4*nn[i]) );
                NEIGH ( i, 4*nn[i]+3 ) = dx[3];
                NI(i,nn[i]++) = j;
                V3DIV ( dx, -dx[3], &NEIGH(j,4*nn[j]) );
                NEIGH ( j, 4*nn[j]+3 ) = dx[3];
                NI(j,nn[j]++) = i;
            }
        }
    pote = 0;
    VZERO (np, X);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( tp[i] == tp[j] )
                if ( tp[i] == Zr_ATOM )
                {
                    if (dx[3] >= SQUARE(RCUT_ZrZr_0)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    pote += 2 * exp( A_ZrZr * (B_ZrZr - dx[3]) + cu );
                    cc = exp( C_ZrZr * (D_ZrZr - dx[3]) + cu );
                    X[i] += cc;
                    X[j] += cc;
                }
                else {}
            else
            {
                if (dx[3] >= SQUARE(RCUT_ZrC_0)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                pote += 2 * exp( A_ZrC * (B_ZrC - dx[3]) +
                                 K_ZrC / (dx[3] - RCUT_ZrC) );
                for (sc=k=0; k<nn[i]; k++)
                { /* screening at Zr */
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(dx, &NEIGH(i,4*k)) / dx[3];
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) ) );
                }
                for (k=0; k<nn[j]; k++)
                { /* screening at C */
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(dx, &NEIGH(j,4*k)) / dx[3];
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) ) );
                }
                sc = sqrt(sc) / exp( K_ZrC / (dx[3] - RCUT_ZrC) -
                                     C_ZrC * dx[3] );
                h = exp( C_ZrC * (D_ZrC - dx[3]) +
                         K_ZrC / (dx[3] - RCUT_ZrC) - sc );
                X[i] += h * h;
                X[j] += h * h;
                pote -= LAMBDA * h;
                if (swapped) SWAP(i,j,k);
            }
        }
    for (i=np; i--;) pote -= (X[i]=sqrt(X[i]));
    VZERO (DIMENSION*np, f);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( tp[i] == tp[j] )
                if ( tp[i] == Zr_ATOM )
                {
                    if (dx[3] >= SQUARE(RCUT_ZrZr_0)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    cv = cu / (dx[3] - RCUT_ZrZr);
                    ff = ( 2 * ( A_ZrZr + cv ) *
                           exp( A_ZrZr * (B_ZrZr - dx[3]) + cu ) -
                           ( C_ZrZr + cv ) *
                           exp( C_ZrZr * (D_ZrZr - dx[3]) + cu ) *
                           ( 0.5 / X[i] + 0.5 / X[j] ) ) / dx[3];
                    V3MUL (ff, dx, force);
                    V3AdD (force, &f[DIMENSION*j]);
                    V3SuB (&f[DIMENSION*i], force);
                }
                else {}
            else
            {
                if (dx[3] >= SQUARE(RCUT_ZrC_0)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                V3DIV ( dx, dx[3], xx );
                for (sc=k=0; k<nn[i]; k++)
                { /* screening at Zr */
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(xx, &NEIGH(i,4*k));
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) ) );
                }
                for (k=0; k<nn[j]; k++)
                { /* screening at C */
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(xx, &NEIGH(j,4*k));
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) ) );
                }
                cu = K_ZrC / (dx[3] - RCUT_ZrC);
                cv = cu / (dx[3] - RCUT_ZrC);
                sc = sqrt(sc) / exp( cu - C_ZrC * dx[3] );
                h = exp( C_ZrC * (D_ZrC - dx[3]) + cu - sc );
                cc = h * h / X[i] + h * h / X[j] + LAMBDA * h;
                ff = ( 2 * ( A_ZrC + cv ) *
                       exp( A_ZrC * (B_ZrC - dx[3]) + cu ) -
                       ( C_ZrC + cv ) * cc * ( 1 + sc ) ) / dx[3];
                V3MUL (ff, dx, force);  /* j from i */
                if (sc < EPS) goto checkout;
                for (k=0; k<nn[i]; k++)
                {
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(xx, &NEIGH(i,4*k));
                    if (tt <= 0) continue;
                    ss = cc / sc / 2 *
                        pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) - cu +
                                    C_ZrC * dx[3] ) );
                    ff = ss * 2 * ( C_ZrC + K_ZrC /
                                    (NEIGH(i,4*k+3) - RCUT_ZrC) /
                                    (NEIGH(i,4*k+3) - RCUT_ZrC) );
                    V3MUL (ff, &NEIGH(i,4*k), gorce);  /* k from i */
                    ss *= BETA_ZrC / tt;
                    V3CROSS ( xx, &NEIGH(i,4*k), m );
                    V3CROSS ( m, &NEIGH(i,4*k), wk );
                    V3CROSS ( xx, m, wx );
                    tt = ss / NEIGH(i,4*k+3);
                    V3ADDmuL (tt, wk, gorce);  /* k from i */
                    V3AdD (gorce, &f[DIMENSION*NI(i,k)]);  /* checkout */
                    V3SuB (&f[DIMENSION*i], gorce);
                    tt = ss / dx[3];
                    V3ADDmuL (tt, wx, force);  /* j from i */
                }
                for (k=0; k<nn[j]; k++)
                {
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(xx, &NEIGH(j,4*k));
                    if (tt <= 0) continue;
                    ss = cc / sc / 2 *
                        pow ( tt /  ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) - cu +
                                    C_ZrC * dx[3] ) );
                    ff = ss * 2 * ( C_ZrC + K_ZrC /
                                    (NEIGH(j,4*k+3) - RCUT_ZrC) /
                                    (NEIGH(j,4*k+3) - RCUT_ZrC) );
                    V3MUL (ff, &NEIGH(j,4*k), gorce);  /* k from j */
                    ss *= BETA_CZr / tt;
                    V3CROSS ( &NEIGH(j,4*k), xx, m );
                    V3CROSS ( m, &NEIGH(j,4*k), wk );
                    V3CROSS ( xx, m, wx );
                    tt = ss / NEIGH(j,4*k+3);
                    V3ADDmuL (tt, wk, gorce);  /* k from j */
                    V3AdD (gorce, &f[DIMENSION*NI(j,k)]);  /* checkout */
                    V3SuB (&f[DIMENSION*j], gorce);
                    tt = ss / dx[3];
                    V3ADDmuL (tt, wx, force);  /* j from i */
                }
              checkout:
                V3AdD (force, &f[DIMENSION*j]);
                V3SuB (&f[DIMENSION*i], force);
                if (swapped) SWAP(i,j,k);
            }
        }
    return;
} /* end potential_force() */


/* Fill in motion, pote, f, stress */
void potential_stress()
{
    register int i, j, k;
    int l, swapped;
    register double cc, ff;
    double h, sc, ss, tt, cu, cv;
    double ds[DIMENSION], dx[4], m[DIMENSION], xx[DIMENSION], wk[DIMENSION],
        wx[DIMENSION], force[DIMENSION], gorce[DIMENSION];
    Config_motion_stress (Config_Aapp_to_Alib, motion, stress);
    IZERO (np, nn);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            if (tp[i] == tp[j]) continue;
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( dx[3] < SQUARE(RCUT_ZrC) )
            {
                dx[3] = sqrt(dx[3]);
                if ((nn[i]==MAX_FIRST) || (nn[j]==MAX_FIRST))
                    pe("potential_stress: MAX_FIRST=%d reached.\n",
                       MAX_FIRST);
                V3DIV ( dx, dx[3], &NEIGH(i,4*nn[i]) );
                NEIGH ( i, 4*nn[i]+3 ) = dx[3];
                NI(i,nn[i]++) = j;
                V3DIV ( dx, -dx[3], &NEIGH(j,4*nn[j]) );
                NEIGH ( j, 4*nn[j]+3 ) = dx[3];
                NI(j,nn[j]++) = i;
            }
        }
    pote = 0;
    VZERO (np, X);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( tp[i] == tp[j] )
                if ( tp[i] == Zr_ATOM )
                {
                    if (dx[3] >= SQUARE(RCUT_ZrZr_0)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    pote += 2 * exp( A_ZrZr * (B_ZrZr - dx[3]) + cu );
                    cc = exp( C_ZrZr * (D_ZrZr - dx[3]) + cu );
                    X[i] += cc;
                    X[j] += cc;
                }
                else {}
            else
            {
                if (dx[3] >= SQUARE(RCUT_ZrC_0)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                pote += 2 * exp( A_ZrC * (B_ZrC - dx[3]) +
                                 K_ZrC / (dx[3] - RCUT_ZrC) );
                for (sc=k=0; k<nn[i]; k++)
                { /* screening at Zr */
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(dx, &NEIGH(i,4*k)) / dx[3];
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) ) );
                }
                for (k=0; k<nn[j]; k++)
                { /* screening at C */
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(dx, &NEIGH(j,4*k)) / dx[3];
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) ) );
                }
                sc = sqrt(sc) / exp( K_ZrC / (dx[3] - RCUT_ZrC) -
                                     C_ZrC * dx[3] );
                h = exp( C_ZrC * (D_ZrC - dx[3]) +
                         K_ZrC / (dx[3] - RCUT_ZrC) - sc );
                X[i] += h * h;
                X[j] += h * h;
                pote -= LAMBDA * h;
                if (swapped) SWAP(i,j,k);
            }
        }
    for (i=np; i--;) pote -= (X[i]=sqrt(X[i]));
    VZERO (DIMENSION*np, f);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( tp[i] == tp[j] )
                if ( tp[i] == Zr_ATOM )
                {
                    if (dx[3] >= SQUARE(RCUT_ZrZr_0)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    cv = cu / (dx[3] - RCUT_ZrZr);
                    ff = ( 2 * ( A_ZrZr + cv ) *
                           exp( A_ZrZr * (B_ZrZr - dx[3]) + cu ) -
                           ( C_ZrZr + cv ) *
                           exp( C_ZrZr * (D_ZrZr - dx[3]) + cu ) *
                           ( 0.5 / X[i] + 0.5 / X[j] ) ) / dx[3];
                    V3MUL (ff, dx, force);
                    V3AdD (force, &f[DIMENSION*j]);
                    V3SuB (&f[DIMENSION*i], force);
                    SYMMAT_ACCUMULATE (force, dx, stress);
                }
                else {}
            else
            {
                if (dx[3] >= SQUARE(RCUT_ZrC_0)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                V3DIV ( dx, dx[3], xx );
                for (sc=k=0; k<nn[i]; k++)
                { /* screening at Zr */
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(xx, &NEIGH(i,4*k));
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) ) );
                }
                for (k=0; k<nn[j]; k++)
                { /* screening at C */
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(xx, &NEIGH(j,4*k));
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) ) );
                }
                cu = K_ZrC / (dx[3] - RCUT_ZrC);
                cv = cu / (dx[3] - RCUT_ZrC);
                sc = sqrt(sc) / exp( cu - C_ZrC * dx[3] );
                h = exp( C_ZrC * (D_ZrC - dx[3]) + cu - sc );
                cc = h * h / X[i] + h * h / X[j] + LAMBDA * h;
                ff = ( 2 * ( A_ZrC + cv ) *
                       exp( A_ZrC * (B_ZrC - dx[3]) + cu ) -
                       ( C_ZrC + cv ) * cc * ( 1 + sc ) ) / dx[3];
                V3MUL (ff, dx, force);  /* j from i */
                if (sc < EPS) goto checkout;
                for (k=0; k<nn[i]; k++)
                {
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(xx, &NEIGH(i,4*k));
                    if (tt <= 0) continue;
                    ss = cc / sc / 2 *
                        pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) - cu +
                                    C_ZrC * dx[3] ) );
                    ff = ss * 2 * ( C_ZrC + K_ZrC /
                                    (NEIGH(i,4*k+3) - RCUT_ZrC) /
                                    (NEIGH(i,4*k+3) - RCUT_ZrC) );
                    V3MUL (ff, &NEIGH(i,4*k), gorce);  /* k from i */
                    ss *= BETA_ZrC / tt;
                    V3CROSS ( xx, &NEIGH(i,4*k), m );
                    V3CROSS ( m, &NEIGH(i,4*k), wk );
                    V3CROSS ( xx, m, wx );
                    tt = ss / NEIGH(i,4*k+3);
                    V3ADDmuL (tt, wk, gorce);  /* k from i */
                    V3AdD (gorce, &f[DIMENSION*NI(i,k)]);  /* checkout */
                    V3SuB (&f[DIMENSION*i], gorce);
                    SYMMAT_accumulate
                        (gorce, NEIGH(i,4*k+3), &NEIGH(i,4*k), stress);
                    tt = ss / dx[3];
                    V3ADDmuL (tt, wx, force);  /* j from i */
                }
                for (k=0; k<nn[j]; k++)
                {
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(xx, &NEIGH(j,4*k));
                    if (tt <= 0) continue;
                    ss = cc / sc / 2 *
                        pow ( tt /  ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) - cu +
                                    C_ZrC * dx[3] ) );
                    ff = ss * 2 * ( C_ZrC + K_ZrC /
                                    (NEIGH(j,4*k+3) - RCUT_ZrC) /
                                    (NEIGH(j,4*k+3) - RCUT_ZrC) );
                    V3MUL (ff, &NEIGH(j,4*k), gorce);  /* k from j */
                    ss *= BETA_CZr / tt;
                    V3CROSS ( &NEIGH(j,4*k), xx, m );
                    V3CROSS ( m, &NEIGH(j,4*k), wk );
                    V3CROSS ( xx, m, wx );
                    tt = ss / NEIGH(j,4*k+3);
                    V3ADDmuL (tt, wk, gorce);  /* k from j */
                    V3AdD (gorce, &f[DIMENSION*NI(j,k)]);  /* checkout */
                    V3SuB (&f[DIMENSION*j], gorce);
                    SYMMAT_accumulate
                        (gorce, NEIGH(j,4*k+3), &NEIGH(j,4*k), stress);
                    tt = ss / dx[3];
                    V3ADDmuL (tt, wx, force);  /* j from i */
                }
              checkout:
                V3AdD (force, &f[DIMENSION*j]);
                V3SuB (&f[DIMENSION*i], force);
                SYMMAT_ACCUMULATE (force, dx, stress);
                if (swapped) SWAP(i,j,k);
            }
        }
    SYMMAT_ACCOMPLISH (H, cc, stress);
    return;
} /* end potential_stress() */


/* Fill in motion, pote, f, stress, v, ep, hc */
void potential_current()
{
    register int i, j, k;
    int l, swapped;
    register double cc, ff;
    double h, sc, ss, tt, cu, cv;
    double dd, ee, pi, pj;
    double ds[DIMENSION], dx[4], m[DIMENSION], xx[DIMENSION], wk[DIMENSION],
        wx[DIMENSION], force[DIMENSION], gorce[DIMENSION];
    V3ZERO (hc);
    Config_motion_stress_v_ep (Config_Aapp_to_Alib, motion, stress, v, ep);
    IZERO (np, nn);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            if (tp[i] == tp[j]) continue;
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( dx[3] < SQUARE(RCUT_ZrC) )
            {
                dx[3] = sqrt(dx[3]);
                if ((nn[i]==MAX_FIRST) || (nn[j]==MAX_FIRST))
                    pe("potential_current: MAX_FIRST=%d reached.\n",
                       MAX_FIRST);
                V3DIV ( dx, dx[3], &NEIGH(i,4*nn[i]) );
                NEIGH ( i, 4*nn[i]+3 ) = dx[3];
                NI(i,nn[i]++) = j;
                V3DIV ( dx, -dx[3], &NEIGH(j,4*nn[j]) );
                NEIGH ( j, 4*nn[j]+3 ) = dx[3];
                NI(j,nn[j]++) = i;
            }
        }
    pote = 0;
    VZERO (np, X);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( tp[i] == tp[j] )
                if ( tp[i] == Zr_ATOM )
                {
                    if (dx[3] >= SQUARE(RCUT_ZrZr_0)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    cc = exp( A_ZrZr * (B_ZrZr - dx[3]) + cu );
                    ep[i] += cc;
                    ep[j] += cc;
                    pote += 2 * cc;
                    cc = exp( C_ZrZr * (D_ZrZr - dx[3]) + cu );
                    X[i] += cc;
                    X[j] += cc;
                }
                else {}
            else
            {
                if (dx[3] >= SQUARE(RCUT_ZrC_0)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                cc = exp( A_ZrC * (B_ZrC - dx[3]) +
                          K_ZrC / (dx[3] - RCUT_ZrC) );
                ep[i] += cc;
                ep[j] += cc;
                pote += 2 * cc;
                for (sc=k=0; k<nn[i]; k++)
                { /* screening at Zr */
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(dx, &NEIGH(i,4*k)) / dx[3];
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) ) );
                }
                for (k=0; k<nn[j]; k++)
                { /* screening at C */
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(dx, &NEIGH(j,4*k)) / dx[3];
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) ) );
                }
                sc = sqrt(sc) / exp( K_ZrC / (dx[3] - RCUT_ZrC) -
                                     C_ZrC * dx[3] );
                h = exp( C_ZrC * (D_ZrC - dx[3]) +
                         K_ZrC / (dx[3] - RCUT_ZrC) - sc );
                X[i] += h * h;
                X[j] += h * h;
                cc = LAMBDA * h;
                ep[i] -= cc / 2.;
                ep[j] -= cc / 2.;
                pote -= cc;
                if (swapped) SWAP(i,j,k);
            }
        }
    for (i=np; i--;)
    {
        X[i] = sqrt(X[i]);
        ep[i] -= X[i];
        pote -= X[i];
    }
    VZERO (DIMENSION*np, f);
    for (i=np; i--;)
        for (l=N->idx[2*i]; l<N->idx[2*i+1]; l++)
        {
            j = N->list[l];
            V3SUB (&s[DIMENSION*j], &s[DIMENSION*i], ds);
            V3ImagE (ds);
            V3M3LENGTH2 (ds, H, dx);
            if ( tp[i] == tp[j] )
                if ( tp[i] == Zr_ATOM )
                {
                    if (dx[3] >= SQUARE(RCUT_ZrZr_0)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    cv = cu / (dx[3] - RCUT_ZrZr);
                    cc = ( A_ZrZr + cv ) *
                        exp( A_ZrZr * (B_ZrZr - dx[3]) + cu );
                    dd = ( C_ZrZr + cv ) *
                        exp( C_ZrZr * (D_ZrZr - dx[3]) + cu );
                    ee = 2 * cc - dd * ( 0.5 / X[i] + 0.5 / X[j] );
                    ff = ee / dx[3];
                    pi = ( cc - dd * 0.5 / X[i] ) / ee;
                    pj = 1 - pi;
                    V3MUL (ff, dx, force);
                    V3AdD (force, &f[DIMENSION*j]);
                    V3SuB (&f[DIMENSION*i], force);
                    cc = force[0] * (pj * v[3*i]  + pi * v[3*j]  ) +
                        force[1] * (pj * v[3*i+1] + pi * v[3*j+1]) +
                        force[2] * (pj * v[3*i+2] + pi * v[3*j+2]);
                    V3ADDmuL (cc, dx, hc);
                    SYMMAT_ACCUMULATE (force, dx, stress);
                }
                else {}
            else
            {
                if (dx[3] >= SQUARE(RCUT_ZrC_0)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                V3DIV ( dx, dx[3], xx );
                cu = K_ZrC / (dx[3] - RCUT_ZrC);
                cv = cu / (dx[3] - RCUT_ZrC);
                ff = 2 * ( A_ZrC + cv ) *
                    exp( A_ZrC * (B_ZrC - dx[3]) + cu ) / dx[3];
                V3MUL (ff, dx, force);  /* j from i */
                V3AdD (force, &f[DIMENSION*j]);
                V3SuB (&f[DIMENSION*i], force);
                SYMMAT_ACCUMULATE (force, dx, stress);
                cc = ( force[0] * (v[3*i]   + v[3*j]  ) +
                       force[1] * (v[3*i+1] + v[3*j+1]) +
                       force[2] * (v[3*i+2] + v[3*j+2]) ) / 2.;
                V3ADDmuL (cc, dx, hc);
                for (sc=k=0; k<nn[i]; k++)
                { /* screening at Zr */
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(xx, &NEIGH(i,4*k));
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) ) );
                }
                for (k=0; k<nn[j]; k++)
                { /* screening at C */
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(xx, &NEIGH(j,4*k));
                    if (tt <= 0) continue;
                    sc += pow ( tt / ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) ) );
                }
                sc = sqrt(sc) / exp( cu - C_ZrC * dx[3] );
                h = exp( C_ZrC * (D_ZrC - dx[3]) + cu - sc );
                if ( LAMBDA == 0 )
                {
                    cc = 1. / X[i] + 1. / X[j];
                    pi = 1. / X[i] / cc;
                    pj = 1 - pi;
                    cc *= h * h;
                }
                else
                {
                    cc = h / X[i] + h / X[j] + LAMBDA;
                    pi = ( h / X[i] + LAMBDA / 2. ) / cc;
                    pj = 1 - pi;
                    cc *= h;
                }
                ff =  - ( C_ZrC + cv ) * cc * ( 1 + sc ) / dx[3];
                V3MUL (ff, dx, force);  /* j from i */
                if (sc < EPS) goto checkout;
                for (k=0; k<nn[i]; k++)
                {
                    if (NI(i,k) == j) continue;
                    tt = 1 + V3DOT(xx, &NEIGH(i,4*k));
                    if (tt <= 0) continue;
                    ss = cc / sc / 2 *
                        pow ( tt / ALPHA_ZrC, BETA_ZrC ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(i,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(i,4*k+3) - cu +
                                    C_ZrC * dx[3] ) );
                    ff = ss * 2 * ( C_ZrC + K_ZrC /
                                    (NEIGH(i,4*k+3) - RCUT_ZrC) /
                                    (NEIGH(i,4*k+3) - RCUT_ZrC) );
                    V3MUL (ff, &NEIGH(i,4*k), gorce);  /* k from i */
                    ss *= BETA_ZrC / tt;
                    V3CROSS ( xx, &NEIGH(i,4*k), m );
                    V3CROSS ( m, &NEIGH(i,4*k), wk );
                    V3CROSS ( xx, m, wx );
                    tt = ss / NEIGH(i,4*k+3);
                    V3ADDmuL (tt, wk, gorce);  /* k from i */
                    V3AdD (gorce, &f[DIMENSION*NI(i,k)]);  /* checkout */
                    V3SuB (&f[DIMENSION*i], gorce);
                    SYMMAT_accumulate
                        (gorce, NEIGH(i,4*k+3), &NEIGH(i,4*k), stress);
                    dd = V3DOT ( gorce, &v[DIMENSION*NI(i,k)] );
                    V3ADDmuL (dd*NEIGH(i,4*k+3), &NEIGH(i,4*k), hc);
                    ee = V3DOT ( gorce, &v[DIMENSION*i] );
                    V3ADDmuL ( pj*(ee-dd), dx, hc);
                    tt = ss / dx[3];
                    V3ADDmuL (tt, wx, force);  /* j from i */
                }
                for (k=0; k<nn[j]; k++)
                {
                    if (NI(j,k) == i) continue;
                    tt = 1 - V3DOT(xx, &NEIGH(j,4*k));
                    if (tt <= 0) continue;
                    ss = cc / sc / 2 *
                        pow ( tt /  ALPHA_CZr, BETA_CZr ) *
                        exp ( 2 * ( K_ZrC / (NEIGH(j,4*k+3) - RCUT_ZrC) -
                                    C_ZrC * NEIGH(j,4*k+3) - cu +
                                    C_ZrC * dx[3] ) );
                    ff = ss * 2 * ( C_ZrC + K_ZrC /
                                    (NEIGH(j,4*k+3) - RCUT_ZrC) /
                                    (NEIGH(j,4*k+3) - RCUT_ZrC) );
                    V3MUL (ff, &NEIGH(j,4*k), gorce);  /* k from j */
                    ss *= BETA_CZr / tt;
                    V3CROSS ( &NEIGH(j,4*k), xx, m );
                    V3CROSS ( m, &NEIGH(j,4*k), wk );
                    V3CROSS ( xx, m, wx );
                    tt = ss / NEIGH(j,4*k+3);
                    V3ADDmuL (tt, wk, gorce);  /* k from j */
                    V3AdD (gorce, &f[DIMENSION*NI(j,k)]);  /* checkout */
                    V3SuB (&f[DIMENSION*j], gorce);
                    SYMMAT_accumulate
                        (gorce, NEIGH(j,4*k+3), &NEIGH(j,4*k), stress);
                    dd = V3DOT ( gorce, &v[DIMENSION*NI(j,k)] );
                    V3ADDmuL ( dd*NEIGH(j,4*k+3), &NEIGH(j,4*k), hc );
                    ee = V3DOT ( gorce, &v[DIMENSION*j] );
                    V3ADDmuL ( pi*(dd-ee), dx, hc );
                    tt = ss / dx[3];
                    V3ADDmuL (tt, wx, force);  /* j from i */
                }
              checkout:
                V3AdD (force, &f[DIMENSION*j]);
                V3SuB (&f[DIMENSION*i], force);
                SYMMAT_ACCUMULATE (force, dx, stress);
                cc = force[0] * (pj * v[3*i]  + pi * v[3*j]  ) +
                    force[1] * (pj * v[3*i+1] + pi * v[3*j+1]) +
                    force[2] * (pj * v[3*i+2] + pi * v[3*j+2]);
                V3ADDmuL (cc, dx, hc);
                if (swapped) SWAP(i,j,k);
            }
        }
    SYMMAT_ACCOMPLISH (H, cc, stress);
    for (i=np; i--;) V3ADDmuL (ep[i], &(v[3*i]), hc);
    return;
} /* end potential_current() */