Ring Polymer tangential Force calculation

LAMMPS LAMMPS Beginners
,
1

Hi LAMMPS developers and users,

I am new to LAMMPS and currently working on a bachelor’s thesis involving simulations of active ring polymers. The system consists of a single ring polymer with FENE bonds and LJ pair interactions.

I would like to apply active forces to exactly two monomers located opposite to each other along the ring. The force on bead i (with position r_i) depends on the local tangential direction, defined as

t_i = r_(i+1) - r_i / |r_(i+1) - r_i|

and and involves both t_i and t_(i+1)​, i.e. the force direction depends explicitly on the positions of the two neighboring beads. Indices are cyclic due to the ring topology.

I have looked into fix addforce, fix active, and related documentation, but I do not see a way to define a force that depends on neighbor geometry in this manner.

My question is: is there a way to implement such a force using existing LAMMPS commands (e.g. via computes or variables), or would this require writing a custom fix in C++ (or possibly using fix python)?

Any guidance or pointers would be greatly appreciated.

Best regards

2

I don’t see a way do implement this without C++ programming, but I have been wrong about things like this before. Finding the previous and next atom from variables would require this to be a global operation and thus require a lot of communication and make the whole process non-parallel and very slow when trying to use variables. Instead, since this involves three atoms that are connected with bonds, I would implement this as a custom angle style. Since you have a bead-spring model, it should not create a conflict with an existing angle potential, and you can pick the desired atoms easily in the topology.

If you have some C++ programming skills, or know somebody that has and is willing to collaborate with you, implementing this by using a simple angle potential source file as template and then rip out the entire block for computing the force and energy and replacing it with your needed force computation. You can keep the loop over the angles and see how to access the three atom positions.