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!
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.
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
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.
I guess no traction on this one…
I implemented the feature and uploaded an MR:
Thanks so much, Berk! Are there any notes on the use of a sinusoidal force applied to a selected group?Again, thanks.
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.
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?
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.
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?
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
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);
real pull_potential(struct pull_t* pull,
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.
Did this topic overstay its welcome? :)
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?
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. :)
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.
That sounds encouraging, I will try messing with it – thanks!
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!
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.
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.