bool NOC;
TermString potential_name;

#define NI(i,j)     ni[(i)*MAX_FIRST+(j)]
#define NEIGH(i,j)  neigh[(i)*MAX_FIRST4+(j)]
static int *nn=NULL, *ni=NULL;
static double *X=NULL, *neigh=NULL;

void potential_realloc (int np)
{
    REALLOC ( potential, nn, np * (1 + MAX_FIRST), int );
    ni = nn + np;
    REALLOC ( potential, X, np * (1 + MAX_FIRST4), double );
    neigh = X + np;
    return;
} /* end potential_realloc() */

double potential (int np, M3 H, double *s, Tp *tp, Neighborlist *N,
                  double *f, M3 stress)
{
    register int i, j, k;
    int l, swapped;
    register double cc, ff;
    double energy, h, sc, ss, tt, cu, cv;
    double ds[DIMENSION], dx[4], m[DIMENSION], xx[DIMENSION], wk[DIMENSION],
        wx[DIMENSION], force[DIMENSION], gorce[DIMENSION];
    potential_realloc(np);
    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: 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;
            }
        }
    
    energy = 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)) continue;
                    dx[3] = sqrt(dx[3]);
                    cu = K_ZrZr / (dx[3] - RCUT_ZrZr);
                    energy += 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)) continue;
                dx[3] = sqrt(dx[3]);
                swapped = (tp[i] == C_ATOM) && (tp[j] == Zr_ATOM);
                if (swapped)
                {
                    SWAP (i,j,k);
                    V3NEG (dx,dx);
                }
                energy += 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;
                energy -= LAMBDA * h;
                if (swapped) SWAP(i,j,k);
            }
        }

    for (i=np; i--;) energy -= (X[i]=sqrt(X[i]));

    VZERO (DIMENSION*np, f);
    SYMMAT_ZERO (stress);
    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)) 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)) 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 (energy);
} /* end potential() */