# Oscillating external force

**URL:** <https://gromacs.bioexcel.eu/t/oscillating-external-force/4000>\
**Category:** User discussions\
**Created:** [April 12, 2022, 7:45am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000 "2022-04-12T07:45:06Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 12, 2022, 7:45am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/1 "2022-04-12T07:45:06Z")

</div>

I think I asked about this before… Is there any obscure way to impose a sinusoidal force acting on a given group of particles in the system? Not through electric fields acting on charged stuff, but as mechanical excitation.

If not, how likely for such a feature request to be successful? :)

Thanks!

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [April 14, 2022, 5:59pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/2 "2022-04-14T17:59:36Z")

</div>

We already had an internal suggestion to implement this. As in the 2022 release we have a mathematical expression parser for pull coordinates, it would be rather easy to add in time as a variable and then the user could do anything time dependent they like.

If you know a bit about coding, it wouldn’t be too hard to add it yourself. We would accept such a feature.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 15, 2022, 9:17am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/3 "2022-04-15T09:17:50Z")

</div>

I have zero knowledge of the Gromacs code, unfortunately. I wish I had the bandwidth to learn it and contribute…  
I was asking if the team could implement it if requested, and you said that this was already discussed internally. If this was released as an early-version undocumented feature, I would be very happy to test it on an actual system. We have an urgent need for it and for now the only realistic option is to recreate the entire system in the uber-slow LAMMPS to do the wiggling…

If someone agreed to do this privately and it yielded the expected phenomena, it would be highly publishable and I am happy to offer coauthorship. Even in the diffusive regime (the least fun, in my opinion), things have already been shown to be rather interesting: [Transport and dispersion across wiggling nanopores | Nature Physics](https://www.nature.com/articles/s41567-018-0239-0)  
I think a long time ago Justin mentioned he might be up for it, but I may remember incorrectly. In any case, it would be a very useful feature.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 19, 2022, 8:43am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/4 "2022-04-19T08:43:20Z")

</div>

I guess no traction on this one…

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [April 20, 2022, 4:22pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/5 "2022-04-20T16:22:41Z")

</div>

I implemented the feature and uploaded an MR:

> **[Add time as a transformation pull coord variable (!2661) · Merge requests ·...](https://gitlab.com/gromacs/gromacs/-/merge_requests/2661)**
>
> This enables arbitrary time dependence for pull coordinates. Time is not allowed as a variable when AWH is active, as AWH assumes the system is in equilibrium. This...

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 20, 2022, 9:52pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/6 "2022-04-20T21:52:29Z")

</div>

Thanks so much, Berk! Are there any notes on the use of a sinusoidal force applied to a selected group?Again, thanks.

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [April 21, 2022, 8:47am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/7 "2022-04-21T08:47:04Z")

</div>

I realized now that you might want to have the force change directly in a sinusoidal fashion. The feature I implemented allow one to change the reference location in any time dependent fashion, but not directly the force. If you choose a relatively small force constant you can get a sinusoidal force if the reference location moves much more than the group you are pulling on.

But I don’t know what you want to achieve.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 21, 2022, 7:43pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/8 "2022-04-21T19:43:27Z")

</div>

It seems like what you have is an umbrella-type pull, in which the ref position is a sinusoid. Indeed, one could achieve oscillations that way, except maybe a bit awkward for my need.  
The simplest way to put what I want would be the constant-force pull applied to a given group, _except_ the force is a sinusoidal function (with amplitude and frequency set exactly), not a constant. Does this make sense?

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [April 21, 2022, 8:23pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/9 "2022-04-21T20:23:36Z")

</div>

Yes, that is what I realized you probably wanted after implementing the feature (which will anyhow be useful for many).

As I said, if your group doesn’t move too much, you could get close to what you want by moving the reference position a lot.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 21, 2022, 8:54pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/10 "2022-04-21T20:54:53Z")

</div>

The group thermally moves in the amounts that are comparable to the displacements I would like to induce externally. :)  
Given that you already have the whole infrastructure in place, would there be any chance implementing the exact thing I need?

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [April 26, 2022, 8:01am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/11 "2022-04-26T08:01:34Z")

</div>

Here is a quick hack on top of release-2022. I hope you can read the diff and apply it manually.

ommit fb4e528e729017ffe5a9f7653a8ed97b7e422bbe (HEAD → pull-oscillatory-hack, release-2022)  
Author: Berk Hess [hess@kth.se](mailto:hess@kth.se)  
Date: Tue Apr 26 09:58:50 2022 +0200

```
Hack for oscillatory COM pulling

Dirty hack to change the potential for constant force pulling
to oscillate in time, using a hard coded period.

```

diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp  
index 5de323006f…5100fbf850 100644  
— a/src/gromacs/pulling/pull.cpp  
+++ b/src/gromacs/pulling/pull.cpp  
@@ -1240,6 +1240,7 @@ static void add\_virial\_coord(tensor vir, const pull\_coord\_work\_t& pcrd, const Pu  
static void calc\_pull\_coord\_scalar\_force\_and\_potential(pull\_coord\_work\_t\* pcrd,  
double dev,  
real lambda,

- 

```
                                                  const double t,
                                                  real* V,
                                                  real* dVdl)

```

{  
@@ -1268,6 +1269,12 @@ static void calc\_pull\_coord\_scalar\_force\_and\_potential(pull\_coord\_work\_t\* pcrd,  
\*dVdl += 0.5 \* dkdl \* gmx::square(dev);  
break;  
case PullingAlgorithm::ConstantForce:  
+#warning “hack to make constant force pulling oscillatory in time”

- 

```
       {

```

- 

```
           const double period = 10.0;

```

- 

```
           const double factor = 1.0 / (2.0 * M_PI * period);

```

- 

```
           k *= std::sin(factor * t);

```

- 

```
       }
       pcrd->scalarForce += -k;
       *V += k * pcrd->spatialData.value;
       *dVdl += dkdl * pcrd->spatialData.value;

```

@@ -1572,7 +1579,7 @@ static void do\_pull\_pot\_coord(const pull\_t& pull,

```
 double dev = get_pull_coord_deviation(pull, pcrd, pbc, t);

```

- calc\_pull\_coord\_scalar\_force\_and\_potential(pcrd, dev, lambda, V, dVdl);

- calc\_pull\_coord\_scalar\_force\_and\_potential(pcrd, dev, lambda, t, V, dVdl);  
}

real pull\_potential(struct pull\_t\* pull,

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [April 26, 2022, 8:27pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/12 "2022-04-26T20:27:28Z")

</div>

Hi Berk,

Thanks a whole lot for keeping me in mind. If period is the actual period, shouldn’t it be  
`const double factor = 2.0 * M_PI / period;` ?

I thought to ask you for the pull.cpp, but the hardcoded period would make using this hack an absolute nightmare, because frequency sweeps are very high on my agenda.

When using constant force,` pull-coord1-rate` seems unused. Can it be used to load the period in this hack?

Again, thank you.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [May 3, 2022, 8:45am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/13 "2022-05-03T08:45:08Z")

</div>

Did this topic overstay its welcome? :)

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [May 3, 2022, 10:08am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/14 "2022-05-03T10:08:27Z")

</div>

Of course it should be 2pi/period (I probably coded this too late at night).

You could use userreal1 from t\_inputrec. Then you would need to add a variable period to t\_pull and set it in init\_pull() from ir. Do you think you can code this yourself?

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [May 3, 2022, 10:31am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/15 "2022-05-03T10:31:40Z")

</div>

I do know C/C++ syntax, if that’s the question. :) I have zero knowledge of this code and with my cowboy mode of hacking my biggest issue here is the fear of introducing something completely unintended on top of the one thing I want.  
On the other hand, I realize that everyone is busy and I certainly don’t feel comfortable asking you to do it. That’s why I hoped the dev team would consider it a welcome feature for _everyone_. Initially, you gave me hope, but it appears that the oscillatory thing you guys had in mind was something else… Alternatively, I am very happy to offer coauthorship in case publishable data results. :)

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [May 3, 2022, 11:11am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/16 "2022-05-03T11:11:28Z")

</div>

It is as simple as adding in struct pull\_t in pull.h:  
real period;  
The in init\_pull() in pull.cpp adding:  
pull-\>period = ir-\>userreal1  
And then using pull.period instead of period in the modification I posted above.  
Then you can use the userreal1 mdp option to set the period.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [May 3, 2022, 7:58pm UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/17 "2022-05-03T19:58:05Z")

</div>

That sounds encouraging, I will try messing with it – thanks!

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [May 4, 2022, 5:40am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/18 "2022-05-04T05:40:02Z")

</div>

I think after some pain I partially succeeded, except pull\_t is in pull\_internal.h, not in pull.h… I added the period as a struct member and the initialization in init\_pull() is also there now.  
Still unclear: how exactly do I access the period inside calc\_pull\_coord\_scalar\_force\_and\_potential()? I couldn’t see any variables pointing to anything of pull\_t type…  
The period is currently hardcoded and that monster seems operational. Thanks!

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [May 4, 2022, 7:42am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/19 "2022-05-04T07:42:06Z")

</div>

Ah, I think I looked at the wrong function then.

I think the easiest solution is to add the period to pull\_coord\_work\_t instead of to t\_pull.

---

<div class="post-metadata">

**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\
**Post date:** [May 4, 2022, 8:12am UTC](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000/20 "2022-05-04T08:12:17Z")

</div>

I won’t be able to do this (new struct members, proper initialization, etc) without outside help, that much is clear to me after looking at the code structure.  
If you could think of a simple way of passing the period to the hack, that will be very much appreciated.

By the way, is there a simple way to only rebuild the one affected executable and not all of gmx if I keep messing with the code? I have zero cmake knowledge/experience, sadly.

[Next page](https://gromacs.bioexcel.eu/t/oscillating-external-force/4000.md?page=2)
