(no subject)

Hello all,

I’m writing to ask for some advice in making a few modifications to the LAMMPS source code, most likely in the form of 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 subsystem, 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 and 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 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 repeated “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 you had a suggestion for how I could address this issue. I realize 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 this 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