Update forces on particles with run-callback

[Daniel Kappe]

Hi all,

I am currently trying to influence magnetic particles in the simulation with a inhomogeneous magnetic field. In that case, the force vector on each particle is dependent on the particles current position. The only way I found to implement this would be a callback function in the run command updating all force vectors every couple of timesteps. But at the moment this does not work and I am not sure if I got the commands right.

This is the update function I came up with, at the moment it just updates a simple position-dependent force

def force_update(timestep):

    for particle in hoomd.group.all():
        force = (-particle.position[1], particle.position[0], 0.0)
        hoomd.data.force_data_proxy(force, particle)
    return 1

And this is my simple initialization. I used the langevin integrator, because I would like to see the effects of brownian motion in my simulations later.

hoomd.context.initialize("")
system = hoomd.init.create_lattice(unitcell=hoomd.lattice.sq(a=1.4), n=15)
all = hoomd.group.all()
N = len(all)

hoomd.md.integrate.mode_standard(dt=0.002)
hoomd.md.integrate.langevin(group=all, kT=0.0, seed=123, noiseless_r=True) 
hoomd.run(2000, callback=force_update, callback_period=10)

 
I would appreciate any kind of advise, especially on how to make it work or how to solve it more elegantly.

Best,
Daniel

[Jens Glaser]

Hi Daniel,

setting forces from python is currently not supported, except if the force is the same on all particles.

See force.constant(). I am working on extending that capability to set the force per particle.

Expect that to be released with HOOMD 2.2

- Jens

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
To post to this group, send email to hoomd...@googlegroups.com.
Visit this group at https://groups.google.com/group/hoomd-users.
For more options, visit https://groups.google.com/d/optout.

[Daniel Kappe]

Hi Jens,

thank you for the quick response. Is there a lower level way to implement force fields, or should I better wait for the release?

Best,
Daniel

[Jens Glaser]

Yeah, sure. You can always implement the functionality on the C++ level, e.g. as a plug-in.

http://hoomd-blue.readthedocs.io/en/stable/developer.html#plugins-and-components

Or if it is generic enough for a wide user base, submit a pull request for a new feature for HOOMD-blue.

Jens

[Joshua Anderson]

Specifically: Since the force is position dependent, write an external potential evaluator for your force. Then you even get GPU and MPI accelerated implementations for free. All you need to write is the form of the function and provide some linking code. See the existing external potentials for examples.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
http://www-personal.umich.edu/~joaander/

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users+unsubscribe@googlegroups.com.
To post to this group, send email to hoomd-users@googlegroups.com.

Visit this group at https://groups.google.com/group/hoomd-users.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users+unsubscribe@googlegroups.com.

[Daniel Kappe]

Hi,

I have now tried to write an external potential evaluator for force fields, but I have problems passing the field to the evaluator. When running a simulation I get a segmentation fault (exit code 139).

I have written a ForceField class, which is bound to python and looks like this

class ForceField {
public:
    void setSize(/*const input of all public parameters*/) {}; // resizes ffield with specified parameters
    void setValue (const unsigned long &t, const unsigned long &m, const unsigned long &n,
                   const unsigned long &l, const Scalar3 &value); // set Values in ffield
    //gets Values from field
    Scalar3 getValue (const unsigned long &t, const unsigned long &m, const unsigned long &n, const unsigned long &l);
    unsigned long _m,_n,_l,_t; // Size of the force field matrix mxnxlxt
    double _ml,_nl,_ll,_tl; // Length after which the force field repeats itself, time interval of repeat

private:
    GPUVector< Scalar3 > ffield;
};

The class functions setSize(), setValue() and getValue are bound to python and can be read and written in python code. For passing to my evaluator, I added a forcefield class to the external.py file, which is a slightly modified version of periodic and e_field

class force_field(_external_force):
    def __init__(self, forcefield ,name=""):
        hoomd.util.print_status_line();

        # initialize the base class
        _external_force.__init__(self, name);

        # create the c++ mirror class
        if not hoomd.context.exec_conf.isCUDAEnabled():
            self.cpp_force = _md.PotentialExternalForceField(hoomd.context.current.system_definition,self.name);
        else:
            self.cpp_force = _md.PotentialExternalForceFieldGPU(hoomd.context.current.system_definition,self.name);

        hoomd.context.current.system.addCompute(self.cpp_force, self.force_name);

        # setup the coefficient options
        self.required_coeffs = None;

        self.field_coeff = None;

        self.forceField = forcefield;

    def process_coeff(self, coeff):
        pass;

    def process_field_coeff(self,field):
        pass;

    def process_forcefield(self, forcefield):
        return forcefield;

And I added a case to the update_coeffs() function like this

def update_coeffs(self):
    [...]
    if self.forceField:
        forcefield = self.process_forcefield(self.forceField);
        self.cpp_force.setForceField(forcefield);

Which should pass the force field to an added function setForceField() in the PotentialExternal.h header

template<class evaluator>
void PotentialExternal<evaluator>::setForceField(field_type field)
{
    m_field = field;
}

everything is then used by my EvaluatorExternalForceField.h header, which is a modified version of EvaluatorExternalPeriodic.h

class EvaluatorExternalForceField
    {
    public:
        //! type of parameters this external potential accepts
        typedef struct param{} param_type;
        typedef class ForceField field_type;

        //! Constructs the constraint evaluator
        /*! \param X position of particle
            \param box box dimensions
            \param forcefield ForceField object containing force field data
        */
        DEVICE EvaluatorExternalForceField(Scalar3 X, const BoxDim& box, const param_type& params,
                                           const field_type& field ,const field_type& forcefield)
            : m_pos(X),
              m_box(box)
            {
                m_field = forcefield;
            }
        /*
          the needs and set functions are just plain copys from EvaluatorExternalPeriodic.h
        */
        DEVICE void evalForceEnergyAndVirial(Scalar3& F, Scalar& energy, Scalar* virial)
            {
                static long int timestep=0;
                long int m,n,l,t; // as previously introduced, the matrix elements
                F.x = Scalar(0.0);
                F.y = Scalar(0.0);
                F.z = Scalar(0.0);
                energy = Scalar(0.0);
             // evaluate matrix elements
                t = (timestep/m_field._tl)%m_field._t;
                m = lround(fmod((m_pos.x) / m_field._ml, m_field._m));
                if (m<0) {m += m_field._m;}
                n = lround(fmod((m_pos.y)/m_field._nl,m_field._n));
                if (n<0) {n += m_field._n;}
                l = lround(fmod((m_pos.z)/m_field._ll,m_field._l));
                if (l<0) {l += m_field._l;}
                F = m_field.getValue(t,m,n,l); // <- The programm crashes on first call of getValue()
                timestep ++;
            }

        #ifndef NVCC
        //! Get the name of this potential
        /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as
            via analyze.log.
        */
        static std::string getName()
            {
            return std::string("forcefield");
            }
        #endif

    protected:
        Scalar3 m_pos;                //!< particle position
        BoxDim m_box;                 //!< box dimensions
        ForceField m_field;             //!< forcefield data
   };

Calls to ForceField class functions work inside the Evaluator as long as they do not try to access the GPUVector elements. In that case I get a segfault. I am not sure how to fix that problem. I am not that good at programming, especially in c++,  therefore I would appreciate any advise.

Thanks,

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.

[Daniel Kappe]

It seems to have cut my post, here is the rest of the code

      }
      DEVICE void evalForceEnergyAndVirial(Scalar3& F, Scalar& energy, Scalar* virial)
         {

             static long int timestep=0; // should give the current timestep
             long int m,n,l,t; // matrix element numbers

             F.x = Scalar(0.0);
             F.y = Scalar(0.0);
             F.z = Scalar(0.0);
             energy = Scalar(0.0);

             t = (timestep/m_field.getTl())%m_field.getT();

             m = lround(fmod((m_pos.x) / m_field._ml, m_field._m));
             if (m<0) {m += m_field._m;}
             n = lround(fmod((m_pos.y)/m_field._nl,m_field._n));
             if (n<0) {n += m_field._n;}
             l = lround(fmod((m_pos.z)/m_field._ll,m_field._l));
             if (l<0) {l += m_field._l;}

             F = m_field.getValue(t,m,n,l); // <- this throws the error

             timestep ++;
         }

     #ifndef NVCC
     //! Get the name of this potential
     /*! \returns The potential name. Must be short and all lowercase, as this is the name energies will be logged as
         via analyze.log.
     */
     static std::string getName()
         {
         return std::string("forcefield");
         }
     #endif

 protected:
     Scalar3 m_pos;                //!< particle position
     BoxDim m_box;                 //!< box dimensions
     ForceField m_field;             //!< forcefield data
};

Everything runs smoothly until the Evaluators tries to access data in the GPUVector. I think it happens at one of the passing steps, but I am not that experienced writing c++ code and pybinds does not make it easier to understand what actually happens. I would appreciate any advise you could give me.

Thanks a lot,

Daniel

                m_field = forcefield<span style="color: #660;" 

[Michael Howard]

I can't really read any of your code, but it looks like the problem is probably in how you are defining your field vs. the potential. The evaluator for the external potential should be the part doing the operations on the particle position. The field should be something like a vector saying which direction things are pointing, and the parameters should be something that is particle-type specific. The external potentials don't currently have a sense of time, so the only way to get this in would be to update the field. Try taking a look at EvaluatorExternal*.h and EvaluatorWalls.h, which might be a good example for you.

Regards,

Mike

[Daniel Kappe]

The forcefield class is just supposed to be a look-up table for the evaluator, the evaluator itself uses code similar to the other evaluators. I could code the calculation of the position dependent force into the evaluator, but I'd prefer passing a precalculated field, as this would be more flexible.

Should I attach the code to a posting?

Best,

Daniel

[Michael Howard]

If you want to use a lookup table, then you probably need to write a ForceCompute, which is a more base-level class, rather than an external potential evaluator. The evaluator you're writing is only a functor, and so I don't think it's possible for you to have memory management for a large tabulated field within that framework.

Regards,

Mike

[Fabio Grillo]

Dear Jens,

will you implement the capability to set the force per particle any time soon?

Thanks,

Fabio

[Michael Howard]

You are now able to set the force per particle group:

http://hoomd-blue.readthedocs.io/en/stable/module-md-force.html#hoomd.md.force.constant

--Mike

[Jens Glaser]

Hi Fabio,

I’ve completed the PR

https://bitbucket.org/glotzer/hoomd-blue/pull-requests/381/add-possibility-of-callback-in/diff

When the unit tests have run successfully and the branch is merged, it should soon be part of release 2.3.

- Jens

[Fabio Grillo]

Thanks a lot Jens for the prompt response and the effort!

I still have a couple of questions, will you also implemented a callback function for the active force class?

Can I use the new force.constant() before it is implemented in the next release? When do you think it will be realesed? 

Last but not least. Is there an efficient way to impose a spatial dependent translational/rotational diffusio=f(x,y)?

Thanks again, you guys should get a medal for all the effort you put in this.

Fabio

[Jens Glaser]

Hi Fabio,

On Apr 16, 2018, at 3:10 PM, Fabio Grillo <fabio.g...@gmail.com> wrote:

Thanks a lot Jens for the prompt response and the effort!

you’re welcome.

I still have a couple of questions, will you also implemented a callback function for the active force class?

I'll leave that as an exercise for now, since I am not currently using active forces myself. Feel free to put up a feature request on bitbucket!

Can I use the new force.constant() before it is implemented in the next release? When do you think it will be realesed? 

Yes, you can checkout the branch force_callback and compile yourself.

I don’t know when exactly Josh will release 2.3, we are polishing a few things up right now.

Last but not least. Is there an efficient way to impose a spatial dependent translational/rotational diffusio=f(x,y)?

That is a higher level question that very much depends on the integrator you’re using and requires one to think about the physical basis and equations first. Currently, hoomd does not implement this, although an obvious candidate would md.integrate.brownian()/md.integrate.langevin(). I don’t know much how one would correctly implement this, but if you have an idea feel free to communicate or implement.

Thanks again, you guys should get a medal for all the effort you put in this.

Feedback appreciated :)

- Jens

[Fabio Grillo]

Dear Jens thanks again for the super quick reply!

About the implementation of spatial dependent traslational/rotational diff. I forgot to mention that I would like to do it using the langevin integrator. If say gamma/gamma_r can be varied in a callback function one could achieve something of this sort. Yet I am not sure one can vary gamma on the fly. Sorry I am an absolute beginner as far as hoomd-blue is concerned. 

Thanks,

Fabio