development of fix for combined MD/minimization run, e.g., for DRUDE/CORE-SHELL SCF

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS-users and -developers,

sorry in advance for the long mail. I am trying to create a fix that would
enable molecular dynamics runs with intermittent minimization of a subset of
degrees of freedom (DOFs) of a system. I.e., minimize subset of DOFs, get
forces, and then integrate Newton’s equations of motion, and repeat. This could
be useful, e.g., for massless core-shell/drude-particle simulations. This
methodology is less efficient than the already implemented
adiabatic/dual-thermostat models, but tends to be more accurate.

This has brought up before on the list, but I haven’t seen any corresponding
implementation, yet. I could really need this methodology and am willing to
spend time working on the implementation, but I would really appreciate if some
of the experts could maybe give me some pointers…

As far as I can see, LAMMPS has all the basic tools/features needed to do the
subsytem minimization without too much hazzle (or so I thought in the
beginning). So no new methods would be necessary, but one would need to combine
the sub-units in LAMMPS in the right fashion. That’s basically what I tried in
the attached source files. For all what follows I used the LAMMPS-ICMS May23rd
2016. My basic workflow/idea within the fix is as follows:

  • create a group that contains degrees of freedom not included in the fixes
    group (call it FRZGRP)
  • create a setforce fix for FRZGRP (call it SETFORCE)
  • register said fix with min_post_force (need to be friends with modify
    class)
  • create an instance of the Min class in update (as done in regular minimize
    run)
  • use that instance to minimize energy
  • also need some contents of the minimize->setup method, hence be friends
    with Min class

For the last couple of points I tried to use as many high level methods as
possible (that’s why I needed to declare FixSubsysMin a friend class to Min and
Modify).

For now the implementation seems to do on a basic level what it’s supposed to be
doing. Minimize once per step and run one MD steps, and repeat. I attached an
example simulation that contains 256 polarizable SPC water models. The
simulation runs stably using the dual-thermostatting, which I used to
pre-equilibrate the structure in the test run. I tried to use my fix together
with different propagators for the equations of motion.

If I do not constrain the water molecules I always run into problems with
polarization catastrophe. If I use fix rigid/nvt/small to integrate the EOM the
system heats up until another catastrophe. If I use fix nvt integration and
shake the simulation crashes even earlier.

I am quite a bit at a loss right now, and would appreciate if someone more
familiar with LAMMPS could take a look and give me an opinion if this fix could
work and where to go on/improve/fix things.

Thank you very much in advance and best regards,

Frank

min.h.patch (316 Bytes)

modify.h.patch (367 Bytes)

fix_subsys_min.cpp (5.82 KB)

fix_subsys_min.h (1.23 KB)

data.pol (126 KB)

tst.in (1.29 KB)

2

Hi - I’m not clear on what you are trying to do.

You can use LAMMPS as-is to alternate between
minimization and a short run (e.g. 1 step or N steps).
By creating a loop
in the input script which has both commands.
And you can perform the minimization on a subset
of atoms (holding others fixed) by using fix setforce.

How is that different than what you are asking?

Steve

3

Hi Steve,

thanks for the response. Okay, my goal is to run MD for explicitly polarizable systems with Drude oscillators (self-consistently optimized, that is). Furthermore, I’d like to be able to use this in conjunction with other packages, in particular the colvars package for free energy calculations.

You’re of course right, one could alternate in an input script between ‘minimize’ and ‘run’ commands. It seemed to me like the less efficient choice, because in each time step both run and minimize would have to be initialized and finalized. Fixes would have to be ‘unfixed’ as well each timestep, dumps would have to be undumped. It’s not that it’s not possible, but it seems like it would accumulate quite a bit of overhead.

Is there a way to access how many steps elapsed during a minimization in the input script? I haven’t seen it and hence the timestep information would be messy.
Only non-compressed/non-binary output options would be possible if one is not interested in the intermittent minimization (which one usually isn’t), because one would need to undump before a minimization and ‘dump_modify append yes’ afterwards. Strided output would seem to be impossible if the elapsed number of timesteps in the minimization is not accessible.

As a last point I would like to use the FIRE minimizer, but that scrambles the velocities which seems to make running dynamics a no-go.

Let me know if this makes sense (or not). Best regards,

Frank

4

For this scenario, I’d suggest one of 2 options. Both with
the goal of not modifying anything in existing LAMMPS.
Maybe there would be a flag or two that need to be
added to Verlet or the Minimizer to enable them to
be driven with new options.

a) write a new “minrun” command, which alternates
every timestep between a minimize and a run (for 1 step).
There are several other commands that perform
both minimizations and runs - see prd, tad. They can
control how much setup is done before after a run or
minimize, by calling lower level methods in Verlet, Min,
etc. They can reset the timestep after a minimization,
so that it appears time is advancing only due to the dynamics,
They can turn off output during the minimizations.
They could control what fixes are applied during minimization
and runs.

b) write a new “fix minimize” or “fix relax” command. It would
be invoked every N time steps during dynamics at the point in
each timestep where
you want the minimization to occur. It would perform
a minimization, with the fixes it needs (e.g. set force), but
it would similarly reset the time, turn off output, etc.

Steve

5

I thought (that I thought) in a similar direction as option b) with the source files/example that I attached in my initial mail. Basically I was trying to invoke a minimization at every timestep as would be appropriate for my case (could add the other “every N timesteps” option later).
Nevertheless, I am stuck right now. The code seems to work in principle, the minimization is called every time step, acts on the subset of degrees of freedom and then one step in the MD integration is done and repeated.
The problem is that the system heats up and particles fall on top of each other which is of course catastrophic for 1/r potentials. I think it is related to something in the rigid body dynamics that I am possibly screwing up. Alternatively, using SHAKE fix leads to similar fails, only much much faster. Not using rigid bodies/constraints is not an option for me. Please note, that the rigid bodies are fixed during the intermittent minimization and are only moved during the dynamics. The system was pre-equilibrated and should run without any problem.

So if you can spare the time, I’d really appreciate it if you could give the code a quick look and tell me if you find something that is obviously missing/wrong. I left some comments tagged by my initials (FU), also where I am unsure. There is also an input for a test case that fails. I understand that it’s quite a bit to ask, I would be happy with a pointer to what might be needed and I’ll try take care of it myself.

Thanks again,

Frank

6

The minimize in LAMMPS does not enforce rigid
body constraints like fix shake or fix rigid.
If you do what I was suggesting as (b), I don’t
think you would need to modify anything in
main LAMMPS (like your first email mentions),
just add a fix. It would instantiate a Min class, etc.

I suggest you try what you have or (b) w/out SHAKE
and see if you can get it to work.

Steve

7

I realize that the minimize doesn’t enforce rigid bodies. The rigid bodies are actually frozen during minimization (because of the setforce fix) and it seemed related to this at first. But similar problems occur for flexible water models.
I am not aware of a flexible water model with Drude oscillators. If I arbitrarily reduce the polarizability of the model I am working with (to make it run stable for at least a short while) I can run it as a flexible model. It heats up like crazy and at some point disintegrates. Temperature seems to stabilize at around 700K before that (but was set to 300K). Again, this is started from a pre-equilibrated snapshot of a normally stably running model. So this shouldn’t happen. I attached the corresponding input. I still don’t see where this issue might come from. I would like to figure out this issue first before starting to rewrite my fix (which in my opinion is option b), but with modifying a tiny bit of the LAMMPS base)

Concerning the not changing anything in LAMMPS. If I would create a new instance of the Min class (called ‘MinI’ in the following) in my fix and would call its run function, then the following methods would be called in the following order:

MinI->run()

MinI->iterate()

MinI->linemin() (depending on minimizer used)

MinI->energy_force()

As I see it, the problem is in the last one, which relies on the modify instance which contains the list of, e.g., post_force things to do (part of modify instance). Hence, if I don’t either write a new energy_force method or create another instance of the modify class within the LAMMPS namespace I don’t see how this should work without my fix being a friend class of modify (because I need to add the setforce fix).
Question: can I easily “overload” the MinI->energy_force() function with one of my own making (which would then not rely on modify, without having to create a derived class or similar)?

So creating a local instance of Min solves one problem, but not everything. I currently don’t see the ‘proper’ solution for not being a friend of modify class, but would appreciate to be enlightened :wink:

Cheers,

Frank

tst.in (1.1 KB)

data.pol (126 KB)

8

As I see it, the problem is in the last one, which relies on the modify instance which contains the list of, e.g., post_force >things to do (part of modify instance). Hence, if I don’t either write a new energy_force method or create another instance of >the modify class within the LAMMPS namespace I don’t see how this should work without my fix being a friend class of >modify (because I need to add the setforce fix).
Question: can I easily “overload” the MinI->energy_force() function with one of my own making (which would then not rely >on modify, without having to create a derived class or similar)?

Sorry, not following this argument. I don’t think your new fix needs to create a new instance

of Min (I may have said that earlier). Update already has a min and integrate instance.

I believe they are always current with whatever min_style settings you’ve made.

So your fix should be able to just invoke methods in update->min.

Whether it’s update->min or a Min class you instantiate, there is only one Modify

class, and Min uses it. The Modify class has lists of all the fixes that min->energy_force()

will invoke, and your fix could be in it (or not). So I’m not clear on why you are asking about

the friend issue. Regardless, adding a one-liner, like a friend setting to some other

LAMMPS class is not an issue, we can do that if needed.

The bigger issue to me is whether you can run dynamics and pause (every step or less often)

and run a minimization to convergence. This is possible at a high level, like the run every command

does, or like the prd and tad commands do. (this is the (a) solution I proposed). But there a higher level command is managing

the invocation and setup of both the Intergrate and Min class methods. By adding

a fix that does a minimization “inside” of a on-going run, that is doing something different.

I’m not certain it will work as (b). There may be some underlying issues I’m not thinking of.

E.g. things like timestep advancement and output triggering will need to be disabled, else the

outer “run” will get screwed up.

Steve

9

As I see it, the problem is in the last one, which relies on the modify
instance which contains the list of, e.g., post_force >things to do (part of
modify instance). Hence, if I don't either write a new energy_force method
or create another instance of >the modify class within the LAMMPS namespace
I don't see how this should work without my fix being a friend class of
>modify (because I need to add the setforce fix).
Question: can I easily "overload" the MinI->energy_force() function with
one of my own making (which would then not rely >on modify, without having
to create a derived class or similar)?

Sorry, not following this argument. I don't think your new fix needs to
create a new instance
of Min (I may have said that earlier). Update already has a min and
integrate instance.
I believe they are always current with whatever min_style settings you've
made.
So your fix should be able to just invoke methods in update->min.

Whether it's update->min or a Min class you instantiate, there is only one
Modify
class, and Min uses it. The Modify class has lists of all the fixes that
min->energy_force()
will invoke, and your fix could be in it (or not). So I'm not clear on why
you are asking about
the friend issue. Regardless, adding a one-liner, like a friend setting to
some other
LAMMPS class is not an issue, we can do that if needed.

The bigger issue to me is whether you can run dynamics and pause (every step
or less often)
and run a minimization to convergence. This is possible at a high level,
like the run every command
does, or like the prd and tad commands do. (this is the (a) solution I
proposed). But there a higher level command is managing
the invocation and setup of both the Intergrate and Min class methods. By
adding
a fix that does a minimization "inside" of a on-going run, that is doing
something different.
I'm not certain it will work as (b). There may be some underlying issues
I'm not thinking of.
E.g. things like timestep advancement and output triggering will need to be
disabled, else the
outer "run" will get screwed up.

well, i'd like to point out that we already have something like that,
with the charge equilibration fixes of the QEQ package, and since the
targeted use case is similar (self-consistent computation of
polarization via core-shell or drude) what is done in fix qeq and its
derivative might be worth a look. this might limit the scope a bit,
but on the other hand, should make things easier to implement and
maintain. when trying to be too general with an implementation, one
might get lost in the details and too many possible conflicts.

axel.

10

@steve: Sorry for the confusion. I think we’re getting on the same track though.
I do NOT want to create a new instance of the Min class. I want to use the existing things in update. This is what I did in the previously attached source files.
My problems might very likely related to what you mentioned in terms of timestep advancement/output triggering… I tried to circumvent these issues by invoking the minimize->ev_set(timestep before minimization) method after the intermittent minimization was successful, plus resetting other counters as well. It does seem to work in the sense that the output and counting correspond to the MD time steps.

Nevertheless, as mentioned before, there is something wrong and my system heats up and equilibrates to a temperature more than twice as high as in the input, irrespective of which temperature control fix used. I currently seem to be unable to find the reason for this and thought you guys might have a clue/insight to where I am screwing up things that I did not notice so far.

Is there something else that needs resetting, besides the timestep, ev_flags, output->next? Those are the ones that I changing/resetting right now.

@axel: thanks for the pointer. I originally wanted to avoid having to write my own minimization algorithm and make use of what LAMMPS already has to offer. I might give it a look, but I’d still prefer the other option, if I can get it working (bigger choice, better tested).
So what would be the major things that might be conflicting between minimization and MD, that you can think of? Also, in terms of temperature compute/control?

Cheers,

Frank

11

[...]

@axel: thanks for the pointer. I originally wanted to avoid having to write
my own minimization algorithm and make use of what LAMMPS already has to
offer. I might give it a look, but I'd still prefer the other option, if I
can get it working (bigger choice, better tested).
So what would be the major things that might be conflicting between
minimization and MD, that you can think of? Also, in terms of temperature
compute/control?

i think you are overestimating the effort to implement a minimizer
algorithm, especially for a rather well behaved system, and are
underestimating the effort to coax a module that was written under the
assumption to be in full control to do only what you want it to do in
a restricted environment.

in terms of understanding your temperature issue, this is difficult to
discuss on such a high-level. you are also not making your life easy
by using a rather complex test system.
have you tried something simpler, e.g. without fix rigid or fix shake?
especially for fix rigid, you have to keep in mind that the fix
assembles rigid body information in the setup phase and then will do
updates of these bodies during time integration. it will only look at
forces and velocities, but only update positions from internal data.
similarly, fix shake has to do some position update extrapolation, and
that can also create trouble (hence the need to have fix shake defined
before nose-hoover time integration fixes).

axel.

12

[…]

@axel: thanks for the pointer. I originally wanted to avoid having to write
my own minimization algorithm and make use of what LAMMPS already has to
offer. I might give it a look, but I’d still prefer the other option, if I
can get it working (bigger choice, better tested).
So what would be the major things that might be conflicting between
minimization and MD, that you can think of? Also, in terms of temperature
compute/control?

i think you are overestimating the effort to implement a minimizer
algorithm, especially for a rather well behaved system, and are
underestimating the effort to coax a module that was written under the
assumption to be in full control to do only what you want it to do in
a restricted environment.

You might be completely right here, and that is why I wanted to consult with you guys…

in terms of understanding your temperature issue, this is difficult to
discuss on such a high-level. you are also not making your life easy
by using a rather complex test system.
have you tried something simpler, e.g. without fix rigid or fix shake?

My latest test used a fully flexible model, still the same issues with temperature.

I also tested a single, in x-direction moving water molecule with one drude oscillator on the oxygen atom. Here, the temperature seems okay (it’s a single molecule…) and the drude particle always relaxes to the position of the oxygen atom, as it should. I now could try to increase the system size in smaller steps, or use a completely different system. I haven’t done that yet.

I just noticed a similarly weird issue while I was looking for another example. If one removes the calculation of long-range electrostatics from dual-nose-hoover-thermostatted simulation with Drude oscillators (see /examples/USER/drude/swm4-ndp), and runs this simulation, the temperature of the atoms relaxes to about two thirds of the input temperature (with thermostat set to input temperature). Not sure if this related. Might actually bring it up with the maintainer of the package.

Frank

13

This is an interesting idea, to look at the damped dynamics variants for fix qeq in the QEQ package.

A couple of these use the FIRE and a standard damped dynamics minimizer (for charge, not for

coordinates). However, they do it thru the Min class, rather they encode the kernel

of those minimizers directly in the fix.

Then you avoid all the timestepping/output issues. Possibly there would be a way

to encapsulate those kernels in the Min class (for its variants) in a way that

a fix like yours could invoke them, more similar to how fix qeq does.

Steve

14

Sorry for re-opening this old post, and sorry for the unexpectedly long delay in
response to your last mails, Steve and Axel.

For reference, and for other users that might find this helpful, I would like to
at least share with you a prototypical implementation of what we were talking
about, i.e., a self-consistent optimization of particle positions during a
molecular dynamics simulation (‘self-consistent particle field (SCPF)’).

I put the code on my github account here: https://github.com/fuulish/SCPF

The implementation is based heavily on the available Min class. I.e., the code
for the fixes is more or less a copy of the respective min commands, but adapted
to be used as fixes that accept groups of particles to be optimized at each MD
step. Like this, it was quite straightforward to create fixes for all the
minimizers available in LAMMPS. The code is very flexible due to this use of
general minimization techniques.

However, the implementation is quite slow. This is due to the required line
search in each iteration. It might well be that using FIRE or QuickMin would
enable competitive speeds (compared to direct SCF algorithms, see below). I
haven’t played much with the parameters for these optimizers to achieve a
reasonable speed.

For special problems, like polarizability described by Drude particles
or inducible point dipoles, it would be more efficient to repeatedly solve the
actual SCF equations directly. I have implemented that in combination with an
always-stable predictor corrector algorithm, which I would much rather like to
contribute to LAMMPS (in particular because of the ASPC algorithm, not the SCF).
I will bother you guys about that in a separate mail, though.

As the SCPF code is basically a copy of the min class(es) I was wondering if it
would make sense to convert some of the methods in the min class to static
methods, which would be callable from a fix like mine? This would mean rewriting
some of the minimization code, but in the end reduce duplicate code. That would
be one way to continue this work, albeit I currently don’t have the time to do
it on my own. If someone is willing to help, be my guest. The code is still
useful, or at least was to me, for example for validation purposes.

Cheers, and thanks again for your help,

Frank