Implementing a field theory interaction add-on


I am looking to add a dynamic field theory package to LAMMPS that would calculate interactions between different types of particles similar to Coulomb interactions. This add-on would compute forces by interpolating different types of particles onto density grids and compute the forces based off the gradients of those density grids. If you could look over this and provide feedback on the proposed implementation, I would be more than grateful for your help.

Right now, the structure I have mentally is that this functionality would be encapsulated in a KSPACE child that would use the PPPM class to put the different particle types on different density grids (A homopolymers and A molecules on an AB diblock get added to the same rho field). The KSPACE child compute function would then call the FFT3d class from the fft3d_wrap.cpp to do the transforms, followed by the multiplication by the FFT of the gradient of the Gaussian potential (The FFT of Gaussian would computed once during the init() of the KSPACE child) and then an inverse Fourier transform using FFT3d class and multiplied by the appropriate interaction parameter (\chi or \kappa) and the average density of the system.

I am not sure how to properly design it so that if the user specifies multiple of these interaction parameters in the input file ((chiAB, chiAC, chiBC, etc.), the code would perform all the necessary interaction calculations before moving on (so as to reduce the amount of convolutions on densities that are needed to calculate). Also, would these gradient fields that exert forces on the particles be considered a force field as defined in 8.5.1 of the manual, and if so, how should that affect where the fields are stored, or should they be stored under the KSPACE child?

I don’t know whether a ‘fix’ or a User package would be better for this implementation. Any suggestions about which one to use would be appreciated as well as what classes should be used for the intermediate steps.



For those interested, my in-house code does the following:

  1. For each molecule, add it to its respective rho field using a particle to mesh scheme (A homopolymers and A molecules on an AB diblock get added to the same rho field).

  2. FFTW-MPI the densities and multiply them by a FFTW-MPI of a Gaussian potential gradient (perform a convolution).

  3. Inverse the convolution back to real space and multiply it by the appropriate \chi parameter and the average density of the system.

  4. Add this result in the gradient field for each monomer type.

  5. Repeat steps 2-4 for every interaction between types of molecules (chiAB, chiAC, chiBC, etc.)

  6. Multiply the accumulated gradients for each field by the grid weight for each monomer.

this is something where it is probably best, if either stan or steve
or both (cc'd) weigh in on.
i only know about this roughly as much as you do (apart from purely
technical issues).