Oscillating external force

User discussions
21

Hi Sasha, once you create a build directory (called build in the GROMACS example pages) you can go back to it after changing source files. Unless you did make clean in there (or equivalent if you used another build driver, e.g. Ninja), the next build will skip re-building of files that don’t need to be updated.

In addition, using ccache is very helpful to keep a cached copy of all object files. To use it, install it based on what your OS requires (Linux is easiest) and use -DGMX_ENABLE_CCACHE=ON). Most automated build tests rely heavily on this feature to run quickly.

Giacomo

22

Thanks! My question though is this: whatever code changes I introduce only affect mdrun. Is it possible to only keep rebuilding this one executable? In other words, is it possible to have gmx mdrun1, gmx mdrun2, etc in the same build directory?

23

I don’t think you can do that, since mdrun is an option to the gmx program, so it would rather be gmx1 mdrun, gmx2 mdrun…

The simplest & safest (i.e. recommended) solution would be to use separate build directories for each version.

Giacomo

24

Yup, that’s what I’ve been doing so far.

25

Is it possible to have more than one pull section in the mdp file? The reason for this question is this: say, I want to calculate the PMF of an ion while the membrane is being shaken as discussed above. Both require pull sections, one being the hacked constant-force applied to the membrane, while the other is zero-velocity pull of that ion. Any chance this could be done without resorting to violence? Thanks!

26

Yes, you can specify multiple biases, which can have different types.

pull-ncoords = 2
pull-coord1-type = umbrella
(rest of pull-coord1-* options)
pull-coord2-type = constant-force
(rest of pull-coord2-* options)
27

Thanks Justin.
I think it unfortunately creates a bit of a mess, because I have a total of three groups involved in two separate biases. Namely, the constant-force pull contains

pull_group1_name = xxx

while in the umbrella pull I have

pull_ngroups = 2
pull_group1_name = yyy
pull_group2_name = zzz

I understand that statements like for example pull_coord1_vec are specific to each bias, but the ones above are not, which seems like a problem. I can post my entire pull section prepared per your suggestion, if needed…

Thanks!

28

Your full input options would be useful. I don’t see why this can’t be done.

29

Here it is:

; Pull 
pull_ncoords         = 2
pull                 = yes
pull-nstxout         = 1000
pull-nstfout         = 1000

;pull 1: hacked constant-force
pull-coord1-type     = constant-force
pull_group1_name     = SHAKING_ATOM
pull_coord1_start    = yes
pull_coord1_geometry = direction-periodic
pull_ngroups         = 1
pull-pbc-ref-prev-step-com = yes
pull-group1-pbcatom  = 1
pull_coord1_groups   = 0 1
pull_coord1_origin   = 0 0 0 ;default
pull_coord1_vec      = 0 0 1
pull_coord1_k        = 1200.0  ;force in kJ/mol/nm units
;end pull code 1

;pull 2: umbrella
pull-coord2-type        = umbrella
pull_ngroups            = 2
pull_group1_name        = TEST_ION
pull_group2_name        = MEMBRANE
pull_coord2_geometry    = cylinder
pull-cylinder-r         = 0.5 
pull-group2-pbcatom     = 487 
pull_coord2_groups      = 2 1 
pull-coord2-vec         = 0 0 1
pull_coord2_rate        = 0.0
pull_coord2_k           = 2000          ; kJ mol^-1 nm^-2
pull_coord2_start       = yes           ; define initial COM distance > 0
;end pull code 2
30

The pull_ngroups setting is global; it applies to any of the biases being applied. So you just define 3 groups that can be used by any bias, then call them by name in the appropriate sections.

31

Is the combination of setting pull_ngroups = 3 along with three separate pull_group*_name statements a correct understanding of your suggestion? In other words,

; Pull 
pull_ncoords         = 2
pull                 = yes
pull-nstxout         = 1000
pull-nstfout         = 1000
pull_ngroups         = 3

;pull 1: hacked constant-force
pull-coord1-type     = constant-force
pull_group1_name     = SHAKING_ATOM
pull_coord1_start    = yes
pull_coord1_geometry = direction-periodic
pull-pbc-ref-prev-step-com = yes
pull-group1-pbcatom  = 1
pull_coord1_groups   = 0 1
pull_coord1_origin   = 0 0 0 ;default
pull_coord1_vec      = 0 0 1
pull_coord1_k        = 1200.0  ;force in kJ/mol/nm units
;end pull code 1

;pull 2: umbrella
pull-coord2-type        = umbrella
pull_group2_name        = TEST_ION
pull_group3_name        = MEMBRANE
pull_coord2_geometry    = cylinder
pull-cylinder-r         = 0.5 
pull-group2-pbcatom     = 487 
pull_coord2_groups      = 3 2 
pull-coord2-vec         = 0 0 1
pull_coord2_rate        = 0.0
pull_coord2_k           = 2000          ; kJ mol^-1 nm^-2
pull_coord2_start       = yes           ; define initial COM distance > 0
;end pull code 2
32

As I understand it, yes, but I’m very simple and have only ever attempted one bias at a time :)

33

I would have been as simple, trust me – if the suggested sinusoidal force hack didn’t use pull code. ;) In any case, thanks a bunch. I will be trying this monster soon enough to report further.

34

The code version with the hack suggested by Berk does this with our pull input – even if the pull_coord1_k for the shaking atom is set to zero. I thought I was close to making it work, but no.

Program: gmx mdrun, version 2022.1
Source file: src/gromacs/mdlib/sim_util.cpp (line 554)
Function: void checkPotentialEnergyValidity(int64_t, const gmx_enerdata_t&, const t_inputrec&)
MPI rank: 0 (out of 4)

Internal error (bug):
Step 0: The total potential energy is -nan, which is not finite. The LJ and
electrostatic contributions to the energy are 17038.4 and 166975,
respectively. A non-finite potential energy can be caused by overlapping
interactions in bonded interactions or very large or Nan coordinate values.
Usually this is caused by a badly- or non-equilibrated initial configuration,
incorrect interactions or parameters in the topology.

For more information and tips for troubleshooting, please check the GROMACS
website at Common Errors — GROMACS webpage https://www.gromacs.org documentation

35

As soon as in the code above the first bias is removed and ncoords etc are adjusted appropriately, the issue disappears, even with the hacked code… This hack seems to turn things upside down when multiple biases are present.

36

Go figure. It appears that crashing stops if the order of the two biases from above is switched and the hacked one is last… I hope Berk could comment on this behavior, so I would get at least some peace of mind.

This leaves us with the following issue for PMF calculations: the umbrella pullf/pullx data with -pf/-px supplied to mdrun will now produce two non-time columns and I only am interested in the first one. Is there a way to either output only the first column? Alternatively, the gmx wham documentation describes adding the groupsel.dat file, which, to the best of my understanding of that description, in my case would contain a bunch of lines as follows

1 0
1 0
1 0
....

Is this correct? Thank you!

37

I don’t understand how my suggested changes could produce NaN or an order effect. If you want me to check you code, please share it.

I don’t know about the details of gmx wham.

38

Oh, I was asking Justin about wham…

Anyway, I just compared actual usable data from hacked and pristine codes and the physics seems to work – much better than expected, actually. We may have introduced something when we [very lightly] messed with amplitude parsing, which we abandoned anyway. I agree that your hack should have caused nothing.
In any case, I don’t think I should bother you with this now, but maybe I will at some point later. Thanks! ;)

p.s. I really wish proper sinusoidal excitation was implemented either as part of pulling or as a separate feature. Please accept the fact that non-bio people love and use Gromacs. I just got back from a nanofluidics conference and there’s a whole bunch of people who would welcome this feature.