Hi everyone,
(This message is a duplicate because the formatting on my previous post got mangled somehow)
I’m writing to ask for some advice in making a few modifications to the LAMMPS source code,
most likely in the form of new fixes. In particular, I’m using LAMMPS to handle the classical
integration in a new quantum-classical methodology. The model I’m studying involves proton
transfer in a polar solvent, where the proton is treated as the relevant quantum mechanical sub-
system, while the rest of the transfer complex and the solvent is classical.
In this approach, what I need LAMMPS to do is:
1.) Provide the electrostatic potential felt by the proton for a given solvent configuration at two
different locations along its reaction path as output at each time step.
2.) Propagate the solvent under the average electrostatic force of the proton in both of these
configurations.
The first item is relatively straightforward, since I can get the potential energy from pre-existing
fixes, and set the proton at the proper location, then use a “run 0” command to get the output I
need, provided I use repeated “run 1” commands to propagate the system one step at a time. The
second item is where the issue arises, since it seems that MD packages don’t typically calculate
this type of average force.
My current idea to solve this is to define a post-force fix which can be called after the forces are
initially calculated to correct them from a single-configuration force to an average. However, in
poking through the code, I noticed that the Verlet::Setup() method doesn’t call the post-force fixes,
and since I’m stepping through the integration with multiple “run 1” commands, this will get
called each time (assuming I understand the program flow correctly), which would leave my first
half Verlet step with incorrect forces.
I was wondering if anyone had a suggestion for how I could address this issue. I realize that using
“pre no” on top of the run command might help, since the correctly-calculated forces from the
Verlet::Run() method should carry over then, but I’ve found that using this command seems to
change the trajectories I get out, at least for my system. The only other ideas I have are to try and
incorporate the potential calculation into the fix and run it automatically at each time step, to get
around the repeated use of “run 1” commands, or to add a line to call my fix in Verlet::Setup(), but
I’d like to avoid modifying core routines if at all possible.
Thanks for your time,
Thomas Allen
Makri Group
UIUC