LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] [EXTERNAL] Re: rolling friction and cohesion models under development for granular models?
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] [EXTERNAL] Re: rolling friction and cohesion models under development for granular models?


From: Eric Murphy <murphyericjames@...24...>
Date: Wed, 15 Nov 2017 13:40:04 -0600

Here's a cohesive model that I used some time ago to simulate rough (non-compliant) micron sized particles. It was adapted from pair_colloid and admittedly is in need of some clean up. This would essentially be the Derjaguin, Mueller, Toporov (DMT) model seen often in literature. Its also used by others in the chemical engineering literature (see refs in paper below). The model does not include granular interactions but can be used in tandem with them via pair_style hybrid.  The input is as follows (resembling pair_style colloid)

dmt Hamaker_constant interatomic_distance contact_radius1 contact_radius2 outer_cutoff

Note that the contact radius may be different from the diameter if you are looking to incorporate surface roughness of a particle.

There are some numerical details about proper time-steps to use for energy conservation in collisions in the paper when combining with Hertz and Hooke models in LAMMPS.  Murphy, E. and Subramaniam, S.  "Binary collision outcomes for inelastic soft-sphere models with cohesion." Powder Tech. 305 (2017) 462-476. (http://www.sciencedirect.com/science/article/pii/S0032591016305769)  Also some hints if you'd like to speed things up, depending on your flow set-up. 

Thought this could be of use for folks to use in the short term, while waiting on the Sandia folks.  I'm also open to folks cleaning it up and incorporating it into LAMMPS if that is desired - meant to do it a long time ago myself.  

HTH.


Regards,
Eric Murphy
Graduate Research Assistant, Mechanical Engineering
Iowa State University


On Mon, Nov 13, 2017 at 2:47 PM, Bruce Ferdowsi <bferdowsi@...3963...> wrote:
Hi, Dan,

Thank you for your reply and information! Great news! I am happy to help with testing and evaluations, if it helps/is needed.

Best wishes,
Bruce


On Nov 13, 2017, at 3:35 PM, Bolintineanu, Dan Stefan <dsbolin@...3...> wrote:

Bruce,
These features are currently still under development. I expect we’ll have a working version of a more general granular pair style that includes rolling friction and cohesion, which we can add to the general release in a few weeks, maybe a couple of months at the latest.
 
Best,
Dan
 
 
From: Axel Kohlmeyer [mailto:akohlmey@...92......] 
Sent: Monday, November 13, 2017 1:01 PM
To: Bruce Ferdowsi <bferdowsi@...36.....3963...>
Cc: LAMMPS Users Mailing List <lammps-users@...6335.....sourceforge.net>; Bolintineanu, Dan Stefan <dsbolin@...3...>
Subject: [EXTERNAL] Re: [lammps-users] rolling friction and cohesion models under development for granular models?
 
 
 
On Mon, Nov 13, 2017 at 2:58 PM, Bruce Ferdowsi <bferdowsi@...3963...> wrote:
Hi, All,
 
I was looking through the materials of the fifth LAMMPS Workshop and Symposium and came across this presentation by Dan S. Bolintineanu:
 

here it is indicated that complex contact potentials for granular simulations, including rolling friction and cohesion models are being added to new versions of lammps. My question is if there is an under development version of lammps with preliminary versions of this new features for public access (e.g. on GitHub) or this is still work in progress for private/(public) use and is not/(yet) released for public?

 
this question is best directed at the author of the presentation (cc'd).
 
axel.
 
 

Thank you and have a great day.

Best wishes,
Bruce


------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...396...sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users



 
-- 
Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.


------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...6297....sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users


/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author: Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_dmt.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_special.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace MathSpecial;

/* ---------------------------------------------------------------------- */

PairDMT::PairDMT(LAMMPS *lmp) : Pair(lmp)
{
  writedata = 1;
}

/* ---------------------------------------------------------------------- */

PairDMT::~PairDMT()
{
  if (allocated) {
    memory->destroy(setflag);
    memory->destroy(cutsq);

    memory->destroy(form);
    memory->destroy(a12);
    memory->destroy(sigma);
    memory->destroy(d1);
    memory->destroy(d2);
    memory->destroy(a1);
    memory->destroy(a2);
    memory->destroy(dr1);
    memory->destroy(dr2);
    memory->destroy(ar1);
    memory->destroy(ar2);
    memory->destroy(diameter);
    memory->destroy(cut);
    memory->destroy(offset);
    memory->destroy(sigma3);
    memory->destroy(sigma6);
    memory->destroy(lj1);
    memory->destroy(lj2);
    memory->destroy(lj3);
    memory->destroy(lj4);
  }
}

/* ---------------------------------------------------------------------- */

void PairDMT::compute(int eflag, int vflag)
{
  int i,j,ii,jj,inum,jnum,itype,jtype;
  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
  double rsq,r,forcelj,factor_lj;
  double r2inv,r6inv,fR,dUR,dUA,phi;
  double K[9],h[4],g[4],c[4];
  int *ilist,*jlist,*numneigh,**firstneigh;

  evdwl = 0.0;
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

  double **x = atom->x;
  double **f = atom->f;
  int *type = atom->type;
  int nlocal = atom->nlocal;
  double *special_lj = force->special_lj;
  int newton_pair = force->newton_pair;

  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;

  // loop over neighbors of my atoms

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    xtmp = x[i][0];
    ytmp = x[i][1];
    ztmp = x[i][2];
    itype = type[i];
    jlist = firstneigh[i];
    jnum = numneigh[i];

    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      factor_lj = special_lj[sbmask(j)];
      j &= NEIGHMASK;

      delx = xtmp - x[j][0];
      dely = ytmp - x[j][1];
      delz = ztmp - x[j][2];
      rsq = delx*delx + dely*dely + delz*delz;
      jtype = type[j];

      if (rsq >= cutsq[itype][jtype]) continue;

      ////////no more cases////////////////////////////
      //switch (form[itype][jtype]) {
      //case SMALL_SMALL:
        //r2inv = 1.0/rsq;
        //r6inv = r2inv*r2inv*r2inv;
        //forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
        //fpair = factor_lj*forcelj*r2inv;
        //if (eflag) evdwl = r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) -
                     //offset[itype][jtype];
        //break;

      //case SMALL_LARGE:
        //c2 = a2[itype][jtype];
        //K[1] = c2*c2;
        //K[2] = rsq;
        //K[0] = K[1] - rsq;
        //K[4] = rsq*rsq;
        //K[3] = K[1] - K[2];
        //K[3] *= K[3]*K[3];
        //K[6] = K[3]*K[3];
        //fR = sigma3[itype][jtype]*a12[itype][jtype]*c2*K[1]/K[3];
        //fpair = 4.0/15.0*fR*factor_lj *
          //(2.0*(K[1]+K[2]) * (K[1]*(5.0*K[1]+22.0*K[2])+5.0*K[4]) *
           //sigma6[itype][jtype]/K[6]-5.0) / K[0];
        //if (eflag)
          //evdwl = 2.0/9.0*fR *
            //(1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) *
             //sigma6[itype][jtype]/K[6]) - offset[itype][jtype];
        //if (rsq <= K[1])
          //error->one(FLERR,"Overlapping small/large in pair DMT");
        //break;

      //case LARGE_LARGE:

        r = sqrt(rsq);
        c[1] = a1[itype][jtype];
        c[2] = a2[itype][jtype];
        c[3] = ar1[itype][jtype];
        c[4] = ar2[itype][jtype];
        K[0] = c[1]+c[2];
        K[1] = c[3]*c[4];
        K[2] = c[3]+c[4];
        K[3] = K[1]/K[2];
        K[4] = sigma[itype][jtype];
        K[5] = r-K[0]+K[4];
        K[6] = powint(K[5],-2);
        K[7] = 1.0/6.0;
        K[8] = powint(K[4],-1);
        //K[8] = 1.0/(K[5]*K[6]);
        //g[0] = powint(K[3],-7);
        //g[1] = powint(K[4],-7);
        //g[2] = powint(K[5],-7);
        //g[3] = powint(K[6],-7);
        //h[0] = ((K[3]+5.0*K[1])*K[3]+30.0*K[0])*g[0];
    	//h[1] = ((K[4]+5.0*K[1])*K[4]+30.0*K[0])*g[1];
    	//h[2] = ((K[5]+5.0*K[2])*K[5]-30.0*K[0])*g[2];
    	//h[3] = ((K[6]+5.0*K[2])*K[6]-30.0*K[0])*g[3];
    	//g[0] *= 42.0*K[0]/K[3]+6.0*K[1]+K[3];
    	//g[1] *= 42.0*K[0]/K[4]+6.0*K[1]+K[4];
    	//g[2] *= -42.0*K[0]/K[5]+6.0*K[2]+K[5];
    	//g[3] *= -42.0*K[0]/K[6]+6.0*K[2]+K[6];

    	//fR = a12[itype][jtype]*sigma6[itype][jtype]/r/37800.0;
    	if (r >= K[0]){
		fpair = -a12[itype][jtype]*K[3]*K[7]*K[6];
        	phi = fpair*K[5];
    	//dUR = phi/r + 5.0*fR*(g[0]+g[1]-g[2]-g[3]);
    	//dUA = -a12[itype][jtype]/3.0*r*((2.0*K[0]*K[7]+1.0)*K[7] +
                                    //(2.0*K[0]*K[8]-1.0)*K[8]);
        	//fpair = phi * K[6]
			;}
    	else {
        	phi = -a12[itype][jtype]*K[3]*K[7]*K[8];
        	fpair = phi*K[8];
		phi *= (1.0+K[8]*(K[0]-r));}
       if (eflag)
          evdwl += phi;


       //////no more cases///////// 
       //break;
      //}
      //////////////////////////

      //not sure why the evdwl takes the factor_lj
      //if (eflag) evdwl *= factor_lj;
      //////////////////

      f[i][0] += delx*fpair;
      f[i][1] += dely*fpair;
      f[i][2] += delz*fpair;
      if (newton_pair || j < nlocal) {
        f[j][0] -= delx*fpair;
        f[j][1] -= dely*fpair;
        f[j][2] -= delz*fpair;
      }

      if (evflag) ev_tally(i,j,nlocal,newton_pair,
                           evdwl,0.0,fpair,delx,dely,delz);
    }
  }

  if (vflag_fdotr) virial_fdotr_compute();
}

/* ----------------------------------------------------------------------
   allocate all arrays
------------------------------------------------------------------------- */

void PairDMT::allocate()
{
  allocated = 1;
  int n = atom->ntypes;

  memory->create(setflag,n+1,n+1,"pair:setflag");
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      setflag[i][j] = 0;

  memory->create(cutsq,n+1,n+1,"pair:cutsq");

  memory->create(form,n+1,n+1,"pair:form");
  memory->create(a12,n+1,n+1,"pair:a12");
  memory->create(sigma,n+1,n+1,"pair:sigma");
  memory->create(d1,n+1,n+1,"pair:d1");
  memory->create(d2,n+1,n+1,"pair:d2");
  memory->create(a1,n+1,n+1,"pair:a1");
  memory->create(a2,n+1,n+1,"pair:a2");
  memory->create(dr1,n+1,n+1,"pair:dr1");
  memory->create(dr2,n+1,n+1,"pair:dr2");
  memory->create(ar1,n+1,n+1,"pair:ar1");
  memory->create(ar2,n+1,n+1,"pair:ar2");
  memory->create(diameter,n+1,n+1,"pair:diameter");
  memory->create(cut,n+1,n+1,"pair:cut");
  memory->create(offset,n+1,n+1,"pair:offset");
  memory->create(sigma3,n+1,n+1,"pair:sigma3");
  memory->create(sigma6,n+1,n+1,"pair:sigma6");
  memory->create(lj1,n+1,n+1,"pair:lj1");
  memory->create(lj2,n+1,n+1,"pair:lj2");
  memory->create(lj3,n+1,n+1,"pair:lj3");
  memory->create(lj4,n+1,n+1,"pair:lj4");
}

/* ----------------------------------------------------------------------
   global settings
------------------------------------------------------------------------- */

void PairDMT::settings(int narg, char **arg)
{
  if (narg != 1) error->all(FLERR,"Illegal pair_style command");

  cut_global = force->numeric(FLERR,arg[0]);

  // reset cutoffs that have been explicitly set

  if (allocated) {
    int i,j;
    for (i = 1; i <= atom->ntypes; i++)
      for (j = i+1; j <= atom->ntypes; j++)
        if (setflag[i][j]) cut[i][j] = cut_global;
  }
}

/* ----------------------------------------------------------------------
   set coeffs for one or more type pairs
------------------------------------------------------------------------- */

void PairDMT::coeff(int narg, char **arg)
{
  if (narg < 6 || narg > 7)
    error->all(FLERR,"Incorrect args for pair coefficients");
  if (!allocated) allocate();

  int ilo,ihi,jlo,jhi;
  force->bounds(arg[0],atom->ntypes,ilo,ihi);
  force->bounds(arg[1],atom->ntypes,jlo,jhi);

  double a12_one = force->numeric(FLERR,arg[2]);
  double sigma_one = force->numeric(FLERR,arg[3]);
  //these are not diameters they are equivalent "diameters" of characteristic asperities
  double d1_one = force->numeric(FLERR,arg[4]);
  double d2_one = force->numeric(FLERR,arg[5]);
  //these are the diameters
  double *radius = atom->radius;

  double cut_one = cut_global;
  if (narg == 7) cut_one = force->numeric(FLERR,arg[6]);

  if (d1_one < 0.0 || d2_one < 0.0)
    error->all(FLERR,"Invalid d1 or d2 value for pair DMT coeff");

  int count = 0;
  for (int i = ilo; i <= ihi; i++) {
    for (int j = MAX(jlo,i); j <= jhi; j++) {
      a12[i][j] = a12_one;
      sigma[i][j] = sigma_one;
      if (i == j && d1_one != d2_one)
        error->all(FLERR,"Invalid d1 or d2 value for pair DMT coeff");
      dr1[i][j] = d1_one;
      dr2[i][j] = d2_one;
      d1[i][j] = 2.0 * radius[i];
      d2[i][j] = 2.0 * radius[j];
      diameter[i][j] = 0.5*(d1[i][j]+d2[1][2]);
      //diameter[i][j] = 0.5*(radius[i]+radius[j]);
      //diameter[i][j] = 0.5*(d1_one+d2_one);
      cut[i][j] = cut_one;
      setflag[i][j] = 1;
      count++;
    }
  }

  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}

/* ----------------------------------------------------------------------
   init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

double PairDMT::init_one(int i, int j)
{
  if (setflag[i][j] == 0) {
    a12[i][j] = mix_energy(a12[i][i],a12[j][j],sigma[i][i],sigma[j][j]);
    sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
    d1[i][j] = mix_distance(d1[i][i],d1[j][j]);
    d2[i][j] = mix_distance(d2[i][i],d2[j][j]);
    dr1[i][j] = mix_distance(dr1[i][i],dr1[j][j]);
    dr2[i][j] = mix_distance(dr2[i][i],dr2[j][j]);
    diameter[i][j] = 0.5 * (d1[i][j] + d2[i][j]);
    cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
  }

  sigma3[i][j] = sigma[i][j]*sigma[i][j]*sigma[i][j];
  sigma6[i][j] = sigma3[i][j]*sigma3[i][j];

///////////////
  //deleted because we have no more cases
  //if (d1[i][j] == 0.0 && d2[i][j] == 0.0) form[i][j] = SMALL_SMALL;
  //else if (d1[i][j] == 0.0 || d2[i][j] == 0.0) form[i][j] = SMALL_LARGE;
  //else form[i][j] = LARGE_LARGE;
//////////////

  // for SMALL_SMALL, a1/a2 do not need to be set
  // for SMALL_LARGE, a1 does not need to be set, a2 = radius for i,j and j,i
  // for LARGE_LARGE, a1/a2 are radii, swap them for j,i

///////////again deleted because we have no more cases
  //if (form[i][j] == SMALL_LARGE) {
    //if (d1[i][j] > 0.0) a2[i][j] = 0.5*d1[i][j];
    //else a2[i][j] = 0.5*d2[i][j];
    //a2[j][i] = a2[i][j];
  //} else if (form[i][j] == LARGE_LARGE) {
    //a2[j][i] = a1[i][j] = 0.5*d1[i][j];
    //a1[j][i] = a2[i][j] = 0.5*d2[i][j];
  //}
////////////////////////
  a2[j][i] = a1[i][j] = 0.5*d1[i][j];
  a1[j][i] = a2[i][j] = 0.5*d2[i][j];
  ar2[j][i] = ar1[i][j] = 0.5*dr1[i][j];
  ar1[j][i] = ar2[i][j] = 0.5*dr2[i][j];
  

  //form[j][i] = form[i][j];
  a12[j][i] = a12[i][j];
  sigma[j][i] = sigma[i][j];
  sigma3[j][i] = sigma3[i][j];
  sigma6[j][i] = sigma6[i][j];
  diameter[j][i] = diameter[i][j];

  double epsilon = a12[i][j]/144.0;
  lj1[j][i] = lj1[i][j] = 48.0 * epsilon * sigma6[i][j] * sigma6[i][j];
  lj2[j][i] = lj2[i][j] = 24.0 * epsilon * sigma6[i][j];
  lj3[j][i] = lj3[i][j] = 4.0 * epsilon * sigma6[i][j] * sigma6[i][j];
  lj4[j][i] = lj4[i][j] = 4.0 * epsilon * sigma6[i][j];

  offset[j][i] = offset[i][j] = 0.0;
  if (offset_flag) {
    double tmp;
    offset[j][i] = offset[i][j] =
      single(0,0,i,j,cut[i][j]*cut[i][j],0.0,1.0,tmp);
  }

  return cut[i][j];
}

/* ----------------------------------------------------------------------
   proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairDMT::write_restart(FILE *fp)
{
  write_restart_settings(fp);

  int i,j;
  for (i = 1; i <= atom->ntypes; i++)
    for (j = i; j <= atom->ntypes; j++) {
      fwrite(&setflag[i][j],sizeof(int),1,fp);
      if (setflag[i][j]) {
        fwrite(&a12[i][j],sizeof(double),1,fp);
        fwrite(&sigma[i][j],sizeof(double),1,fp);
        fwrite(&dr1[i][j],sizeof(double),1,fp);
        fwrite(&dr2[i][j],sizeof(double),1,fp);
        fwrite(&cut[i][j],sizeof(double),1,fp);
      }
    }
}

/* ----------------------------------------------------------------------
   proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairDMT::read_restart(FILE *fp)
{
  read_restart_settings(fp);
  allocate();

  int i,j;

  for (i = 1; i <= atom->ntypes; i++)
    for (j = i; j <= atom->ntypes; j++) {
      if (comm->me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
      if (setflag[i][j]) {
        if (comm->me == 0) {
          fread(&a12[i][j],sizeof(double),1,fp);
          fread(&sigma[i][j],sizeof(double),1,fp);
          fread(&dr1[i][j],sizeof(double),1,fp);
          fread(&dr2[i][j],sizeof(double),1,fp);
          fread(&cut[i][j],sizeof(double),1,fp);
        }
        MPI_Bcast(&a12[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&dr1[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&dr2[i][j],1,MPI_DOUBLE,0,world);
        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
      }
    }
}

/* ----------------------------------------------------------------------
   proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairDMT::write_restart_settings(FILE *fp)
{
  fwrite(&cut_global,sizeof(double),1,fp);
  fwrite(&offset_flag,sizeof(int),1,fp);
  fwrite(&mix_flag,sizeof(int),1,fp);
}

/* ----------------------------------------------------------------------
   proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairDMT::read_restart_settings(FILE *fp)
{
  int me = comm->me;
  if (me == 0) {
    fread(&cut_global,sizeof(double),1,fp);
    fread(&offset_flag,sizeof(int),1,fp);
    fread(&mix_flag,sizeof(int),1,fp);
  }
  MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}

/* ----------------------------------------------------------------------
   proc 0 writes to data file
------------------------------------------------------------------------- */

void PairDMT::write_data(FILE *fp)
{
  for (int i = 1; i <= atom->ntypes; i++)
    fprintf(fp,"%d %g %g %g %g\n",i,a12[i][i],sigma[i][i],dr1[i][i],dr2[i][i]);
}

/* ----------------------------------------------------------------------
   proc 0 writes all pairs to data file
------------------------------------------------------------------------- */

void PairDMT::write_data_all(FILE *fp)
{
  for (int i = 1; i <= atom->ntypes; i++)
    for (int j = i; j <= atom->ntypes; j++)
      fprintf(fp,"%d %g %g %g %g %g\n",i,
              a12[i][j],sigma[i][j],dr1[i][j],dr2[i][j],cut[i][j]);
}

/* ---------------------------------------------------------------------- */

double PairDMT::single(int i, int j, int itype, int jtype, double rsq,
                           double factor_coul, double factor_lj,
                           double &fforce)
{
  double K[9],h[4],g[4],c[4];
  double r,r2inv,r6inv,forcelj,phi,fR,dUR,dUA;

  ////////////////////no more cases
  //switch (form[itype][jtype]) {
  //case SMALL_SMALL:
    //r2inv = 1.0/rsq;
    //r6inv = r2inv*r2inv*r2inv;
    //forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
    //fforce = factor_lj*forcelj*r2inv;
    //phi = r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) -
      //offset[itype][jtype];
    //break;

  //case SMALL_LARGE:
    //c2 = a2[itype][jtype];
    //K[1] = c2*c2;
    //K[2] = rsq;
    //K[0] = K[1] - rsq;
    //K[4] = rsq*rsq;
    //K[3] = K[1] - K[2];
    //K[3] *= K[3]*K[3];
    //K[6] = K[3]*K[3];
    //fR = sigma3[itype][jtype]*a12[itype][jtype]*c2*K[1]/K[3];
    //fforce = 4.0/15.0*fR*factor_lj *
      //(2.0*(K[1]+K[2])*(K[1]*(5.0*K[1]+22.0*K[2])+5.0*K[4]) *
       //sigma6[itype][jtype]/K[6] - 5.0)/K[0];
    //phi = 2.0/9.0*fR *
      //(1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) *
       //sigma6[itype][jtype]/K[6]) - offset[itype][jtype];
    //break;

  //case LARGE_LARGE:
  /////////////////////////////////////////     
    r = sqrt(rsq);
    c[1] = a1[itype][jtype];
    c[2] = a2[itype][jtype];
    c[3] = ar1[itype][jtype];
    c[4] = ar2[itype][jtype];
    K[0] = c[1]+c[2];
    K[1] = c[3]*c[4];
    K[2] = c[3]+c[4];
    K[3] = K[1]/K[2];
    K[4] = sigma[itype][jtype];
    K[5] = r-K[0]+K[4];
    K[6] = powint(K[5],-2);
    K[7] = 1.0/6.0;
    K[8] = powint(K[4],-1);
    //K[8] = 1.0/(K[5]*K[6]);
    //g[0] = powint(K[3],-7);
    //g[1] = powint(K[4],-7);
    //g[2] = powint(K[5],-7);
    //g[3] = powint(K[6],-7);
    //h[0] = ((K[3]+5.0*K[1])*K[3]+30.0*K[0])*g[0];
    //h[1] = ((K[4]+5.0*K[1])*K[4]+30.0*K[0])*g[1];
    //h[2] = ((K[5]+5.0*K[2])*K[5]-30.0*K[0])*g[2];
    //h[3] = ((K[6]+5.0*K[2])*K[6]-30.0*K[0])*g[3];
    //g[0] *= 42.0*K[0]/K[3]+6.0*K[1]+K[3];
    //g[1] *= 42.0*K[0]/K[4]+6.0*K[1]+K[4];
    //g[2] *= -42.0*K[0]/K[5]+6.0*K[2]+K[5];
    //g[3] *= -42.0*K[0]/K[6]+6.0*K[2]+K[6];

    //fR = a12[itype][jtype]*sigma6[itype][jtype]/r/37800.0;
    if (r >= K[0]){
	fforce = -a12[itype][jtype]*K[3]*K[7]*K[6];
        phi = fforce*K[5];
    //dUR = phi/r + 5.0*fR*(g[0]+g[1]-g[2]-g[3]);
    //dUA = -a12[itype][jtype]/3.0*r*((2.0*K[0]*K[7]+1.0)*K[7] +
                                    //(2.0*K[0]*K[8]-1.0)*K[8]);
        ;}
    else {
        phi = -a12[itype][jtype]*K[3]*K[7]*K[8];
        fforce = phi*K[8];
	phi *= (1.0+K[8]*(K[0]-r));}

    //phi += a12[itype][jtype]/6.0*(2.0*K[0]*(K[7]+K[8])-log(K[8]/K[7])) -
      //offset[itype][jtype];
    //break;
  //}

  return phi;
}
/* -*- c++ -*- ----------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#ifdef PAIR_CLASS

PairStyle(dmt,PairDMT)

#else

#ifndef LMP_PAIR_DMT_H
#define LMP_PAIR_DMT_H

#include "pair.h"

namespace LAMMPS_NS {

class PairDMT : public Pair {
 public:
  PairDMT(class LAMMPS *);
  virtual ~PairDMT();
  virtual void compute(int, int);
  void settings(int, char **);
  void coeff(int, char **);
  double init_one(int, int);
  void write_restart(FILE *);
  void read_restart(FILE *);
  void write_restart_settings(FILE *);
  void read_restart_settings(FILE *);
  void write_data(FILE *);
  void write_data_all(FILE *);
  double single(int, int, int, int, double, double, double, double &);

 protected:
  enum {SMALL_SMALL,SMALL_LARGE,LARGE_LARGE};

  double cut_global;
  double **cut;
  double **a12,**d1,**d2,**diameter,**a1,**a2,**offset,**ar1,**ar2,**dr1,**dr2;
  double **sigma,**sigma3,**sigma6;
  double **lj1,**lj2,**lj3,**lj4;
  int **form;

  void allocate();
};

}

#endif
#endif

/* ERROR/WARNING messages:

E: Overlapping small/large in pair colloid

This potential is infinite when there is an overlap.

E: Overlapping large/large in pair colloid

This potential is infinite when there is an overlap.

E: Illegal ... command

Self-explanatory.  Check the input script syntax and compare to the
documentation for the command.  You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.

E: Incorrect args for pair coefficients

Self-explanatory.  Check the input script or data file.

E: Invalid d1 or d2 value for pair colloid coeff

Neither d1 or d2 can be < 0.

*/