Invoking force_clear() from a fix?

LAMMPS LAMMPS Mailing List Mirror
1

Hello, I am writing a custom fix and trying to invoke force_clear(), located in verlet.cpp.

I first tried the following line, which did not work:

update->integrate->verlet->force_clear()

After some searching, I found many instances of the following line. I assume that some kind of “find” based on the integrate_style name would need to be applied to the earlier line, but am not sure of the exact syntax.

if (strstr(update->integrate_style, “verlet”))

I also tried including verlet.h in the fix, and simply invoking force_clear(), but that also did not work.

Any help is appreciated.

Thank you,
Anne

2

The force_clear member in verlet is protected so you wouldn’t be able to invoke it directly. You can try making it public for your purpose or just copy paste the functionality of force_clear() as a member function in your fix.

Sincerely,
Adrian Diaz

3

Thank you Adrian,

I tried to do the last option (copy-pasting the function into the fix), but it is leading to undefined errors and segmentation faults. Since everything is already built into the verlet structure, I was thinking it would be more straightforward to make a spinoff of verlet (verlet/custom), including a public force_clear_custom() function (basically a duplicate of the protected force_clear() )

From a fix, would one invoke the public member function of a run_style? I have some preliminary syntax below, but am not quite sure how to put it together.

if (strstr(update->integrate_style, “verlet”))
update->integrate->?->force_clear_custom()

4

How much of the functionality of the force clear method are you looking to replicate? and why? Making a whole copy of verlet seems like a very roundabout approach to me. Note that most of the variables and functions accessed and invoked in the force clear method don’t even belong the verlet class. You just need to obtain the values of those flags being checked if your application is concerned with those. As for the segfaults, I’m not sure where in your fix this is being invoked and whether the force clear method would even be the cause, it may be due to having all the flags uninitialized in your fix if you literally copy pasted the code or just trying to memset arrays that haven’t been allocated yet.

5

I think you need to refresh you C++ knowledge.

Hello, I am writing a custom fix and trying to invoke force_clear(), located in verlet.cpp.

I first tried the following line, which did not work:

update->integrate->verlet->force_clear()

This cannot work. the Verlet class is derived from Integrate it. The integrate class does not have a member called “verlet”.
This is standard C++ polymorphism. Actual Integrator classes are derived from the Integrator class which defines the API and through several “pure” functions (i.e. virtual functions initialized to a NULL pointer) determines which methods must be implemented in the derived class. the upper level class, Update, is a composite that includes an “integrate” and a “minimize” class pointer, which then have to be populated by one of the available classes derived from Integrate or Minimize, respectively. The choice is done through the “run_style” and “min_style” commands.

After some searching, I found many instances of the following line. I assume that some kind of “find” based on the integrate_style name would need to be applied to the earlier line, but am not sure of the exact syntax.

if (strstr(update->integrate_style, “verlet”))

if you do something like that, and you will also have to add a test, that your code will only work with run style verlet.
also, this test should no longer be done this way, but you should include the utils.h header and use
if (utils::strmatch(update->integrate_style,"^verlet")) { …

I also tried including verlet.h in the fix, and simply invoking force_clear(), but that also did not work.

of course not. same thing. that is not how C++ works.

as adrian already mentioned, you either have to make the Verlet::force_clear() method public, or better in this case, declare your new class as a friend class in the Verlet header.

then you have to do (after checking that you really do have a Verlet integrator class or a class derived from it):

Verlet *verlet = (Verlet *) update->integrate;

and after that you can use:

verlet->force_clear();

whether that approach overall is a good idea, i cannot say, but i doubt it. however with the lack of specific information, i cannot say anything definite.

axel.

6

Hi Axel,

Thank you for the clarification of the class structure/syntax; that helped a lot. I will look into the friend class.

To give context, in the custom fix, I want to compute the total potential energy of 2 states

  • State 1: the existing system during the current timestep
  • State 2: same configuration, but with a different charge distribution calculated by the fix

The fix then decides, based on the energy, with which configuration to proceed. The one part I was stuck on was clearing out the forces, to avoid accumulation over the multiple calculations. The surrounding code is below:

// State 1
[ make and save atomic-vector changes]
[ apply relevant reverse and forward comm to update atomic info on all processors ]

//Independent force computation for state 1, copied from verlet.cpp
force_clear(); // having a problem invoking this

if (force->pair) force->pair->compute(eflag,vflag);
if (force->kspace) force->kspace->compute(eflag,vflag);
if (force->newton) comm->reverse_comm();
// Independent PE computation for state 1, copied from compute_pe.cpp

one = 0.0; scalar1 = 0.0;
if (force->pair) one += force->pair->eng_vdwl + force->pair->eng_coul;
MPI_Allreduce(&one,&scalar1,1,MPI_DOUBLE,MPI_SUM,world);
if (force->kspace) scalar1 += force->kspace->energy;

[ do a similar force/pe computation for state 2]

[make a decision to proceed with state 1 or 2 based on scalar1 vs. scalar 2]

Following up on Adrian’s email, I am leaning more towards invoking the method from verlet than the copy-paste approach. However, this is a summary of how it was copy-pasted and what errors I got, for reference.

  • I deactivated the torque flag check, because I know for sure I have rigid bodies and torques in my system

  • I am not sure what the extraflag means

  • I could not compile atom->avec->force_clear in the fix (got an error).

…/fix_qeq_shielded_ecam.cpp:181:28: error: invalid use of incomplete type ‘class LAMMPS_NS::AtomVec’
if (extraflag) atom->avec->force_clear(0,nbytes);

  • This was the modified copy-pasted code:

size_t nbytes;

//if (external_force_clear) return;

// clear force on all particles
// if either newton flag is set, also include ghosts
// when using threads always clear all forces.

int nlocal = atom->nlocal;

if (neighbor->includegroup == 0) {
nbytes = sizeof(double) * nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;

if (nbytes) {
memset(&atom->f[0][0],0,3nbytes);
//if (torqueflag)
memset(&atom->torque[0][0],0,3nbytes);
if (extraflag) atom->avec->force_clear(0,nbytes);
}

// neighbor includegroup flag is set
// clear force only on initial nfirst particles
// if either newton flag is set, also include ghosts

} else {
nbytes = sizeof(double) * atom->nfirst;

if (nbytes) {
memset(&atom->f[0][0],0,3nbytes);
//if (torqueflag)
memset(&atom->torque[0][0],0,3nbytes);
if (extraflag) atom->avec->force_clear(0,nbytes);
}

if (force->newton) {
nbytes = sizeof(double) * atom->nghost;

if (nbytes) {
memset(&atom->f[nlocal][0],0,3nbytes);
//if (torqueflag)
memset(&atom->torque[nlocal][0],0,3nbytes);
if (extraflag) atom->avec->force_clear(nlocal,nbytes);
}
}

7

Hmmm… what you are doing seems like you are inviting a lot of problems, e.g. from fix ordering and early or late modifications during time integration.

this seems to me much more natural to implement as a multi-partition calculation (like parallel tempering, or parallel replica dynamics).
in fact, i see some similarities to the flow of control in PRD.
to summarize, in PRD you start from the same configuration on multiple partitions, decorrelate those and then run decorrelated simulation waiting for a rare event to occur in one of them. if that happens, you do a check and then pick the configuration where the event happened and continue the procedure from there.
in your case you would do the decision/synchronization in every step, so the control loop could be essentially that of the run command, but just at some point (post_integrate()?) do your modifications, and then at a second point (post_force() or end_of_step()?) decide which variant to follow. your fix would then apply the modifications and check which is the preferred variant and synchronize coordinates and velocities across all partitions. since you seem to be doing this in every step, you would not need a separate command like prd or temper as those need to run independently for a given number of steps. this would easily to quite a few partitions before the need for collective operations will limit scaling.

another, also IMO cleaner approach would be to implement a custom run style instead of a fix, i.e. build your own custom variant of verlet instead of subsuming that functionality into a fix. what you are doing is effectively modifying the verlet loop, but you are trying to do it from a fix, which is bound to lead to all kinds of trouble since you are mixing different “levels” of control.
inside your custom verlet loop you can do something similar as is done with run-style respa. in that, you need to keep copies of forces from the different RESPA levels and then compute the total force at each step by summing over them. you could implement a similar copy for your force/velocity/coordinate data from one possible configuration, then reset and go back to do the next variant and then do your magic that decides which one to take and copy or not and continue with the next loop. this would not require multiple partitions, but you need to keep local copies (but only O(nlocal)). you would have to be careful with other computes or fixes that do analysis since those may need to be run multiple times and that should be avoid, since quite a few assume they are only run once per time step.

[…]

Hi Axel,

[…]

  • I could not compile atom->avec->force_clear in the fix (got an error).

…/fix_qeq_shielded_ecam.cpp:181:28: error: invalid use of incomplete type ‘class LAMMPS_NS::AtomVec’
if (extraflag) atom->avec->force_clear(0,nbytes);

again, this is a trivial issue and confirms that you need to refresh your C++ knowledge. the atom.h header includes a forward declaration of the AtomVec class (the base class of all atom styles). that is sufficient for as long as you only pass a pointer to a class instance, but when you need to access the content of the class you need to include the full class definition from the atom_vec.h header file.

axel.

8

Hi Axel, something like that did come to mind while making the verlet_custom copy, but I was not sure if it would get too hackish. I will look into it as well. Thank you for the thorough feedback on that idea.

9

axel.