New force filed implementaion

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.

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:

If it helps, I thought the beginning of this paper was pretty helpful
for explaining how quaternions represent 3D rotations:

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.