It seems like MC would be the possible way to go ahead at this moment with this potential.

I read a recent paper where they use WCA and Yukawa potentials for janus ellipsoids particles.

http://www.nature.com/nmat/journal/v14/n1/full/nmat4111.html

In this WCA potential the only missing part is the angle dependence of patch interactions (which is in Kern-Frenkel model).

You realize of course, that there's no reason you are forced to use

the Kern-Frenkel model to get angle-dependence.

Actually, it looks like they are not using simple WCA pairwise attraction.

They multiplied (divided) the WCA potential by function which depends

on the angle. (See equation (3) from the supplemental materials of

that paper).

In any case, in this paper they are using a smooth continuous

potential to model the interactions between particles. And they are

using molecular dynamics (langevin/brownian dynamics) to run the

simulation (not Monte Carlo).

But you will have to add a new pair style to LAMMPS.

Take a look at the code in the "PairGayBerne::compute()" function.

(lines ~87-230 from the "pair_gayberne.cpp" located in the

"src/ASPHERE" subdirectory distributed with LAMMPS). If you can

understand the code well enough to modify it to implement the function

you're interested in, then this project will be feasible for you.

You will need to understand what quaternions are before you attempt to

read or edit this code. (You will have to learn that on your own.)

Here's an obligatory link to wikipedia:

http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation

If it helps, I thought the beginning of this paper was pretty helpful

for explaining how quaternions represent 3D rotations:

http://scripts.iucr.org/cgi-bin/paper?S0108767387010535

If you want to simulate a patchy surface with tetragonal symmetry,

...then dream up a function with tetragonal angle dependence, and

multiply it by your radial dependence. Then modify

"pair_gayberne.cpp" so that it calculates this function.

But, I also want to have the interaction between non-patchy particles.

For example, I have a big particle with four patches on it (with a tetrahedral symmetry).

I would like to have interactions between patches-patches, patches-particles and particles-particles.

If you have two different kinds of particles, then you will have to

dream up different functions for the interactions between the various

kinds of particles in your simulation.

Please let me know if it is possible to do

Sure, but not trivial. You will have to edit some LAMMPS pair_style

code that probably does funny stuff with quaternions.

Please let me know if there are other potentials where I can have all the interactions (as I mention above) with the angle dependence on patch interactions.

As for coming up with a good function(s), that's the fun part. You

need to be the one to pick that.

I don't think there's anything special about the function they used in

the Nature-Methods paper. (I'm sure they just made it up, like

everyone else does.)

Again, if you are willing to build the patchy-particles from many

smaller atoms (similar to the models in the pictures I sent, and the

approach Steve talked about), then this should be comparatively easy

to simulate in LAMMPS (without dealing with quaternions or making any

modifications to the LAMMPS source code). You can just use ordinary

Lennard-Jones interactions between these smaller atom particles for

example, and use fix langevin. This is by far the easiest choice.

Cheers

Andrew