USER-COLVAR Package: Extended Lagrangian and Langevin coupling to atomic DOF

LAMMPS LAMMPS Mailing List Mirror
1

Greetings,

I am currently using the USER-COLVAR package to implement a method similar to Temperature Accelerated MD. In the attached simple example, there are two systems, both are identical except one uses the extended Lagrangian with a Langevin thermostat on the extended DOF (the C-C-C-C dihedral angle in n-butane), the control uses the same CV, but the extended Lagrangian section of the colvars input file is commented out. In this way, I can still use the USER-COLVARS package to track the value of this CV.

What I would expect to happen is that the system with the extended Lagrangian with a Langevin thermostat would differ from the control simulation. I.e., the extended DOF would be thermostated at a higher temperature, which would be coupled by a harmonic spring to the atomic DOF, causing the dihedral angle in the system with the extended Lagrangian to sample the dihedral angle differently from the control.

What does happen is that the value for the dihedral angle is identical (to the last decimal place) between these two simulations. In fact, all atomic coordinates and thermo output are identical between the two systems. The value for the extended DOF dihedral angle (in the colvar trajectory file) does differ from the actual system, however, it appears that there is no force fed back to the atomic DOFs from the extended DOFs, to which they should be coupled.

Perhaps, I have made an error in syntax or my understanding of how this functionality works… I greatly appreciate any clarification on why this is happening… is this expected behavior?

I am happy to send the simulation output for anyone that wishes to take a closer look, but it was large (7.3 MB) and so I have only attached the input files to this email.

P.S. I am using ReaxFF, which is strange for this example system, however the actual system where this technique will be used is reactive… this is just a greatly simplified example that I have constructed to illustrate the issue.

Best regards, Spence

lammps_small_example_input.tar.gz (6.08 KB)

2

Hi Spencer, can you check if this problem persists using an alternative test system? It seems like your dihedral is stuck very close to 180°, and I’m wondering whether that is an issue.

I just re-ran a test using a dihedral on a softer model system (specifically, the deca-alanine peptide used in NAMD):
https://github.com/colvars/colvars/tree/master/namd/tests/library/000_dihedral_harmonic-fixed

and I do see different behavior with/without the extended mass.

Giacomo

3

Hi Giacomo,

Thanks for taking a look. I created a new example with the alanine peptide, still using reaxFF. I used minimal output so that I could include both the input and output files in the attachment, this time.

As you can see, the trajectories with/without the extended mass are identical to the last decimal place. If the trajectories were close to identical, but differed by, e.g., 0.1%… then I would suspect that the parameters I was using for the extended system were not optimal. However, since the two systems are exactly identical in every way that I can tell (coordinates, qeq charges, value of the dihedral, energy, pressure, temperature), I am concerned that there is simply no force coupled to the atomic dof from the extended dof.

Not sure if you have seen USER-COLVARS used with reaxff… maybe it is a bug related to this??

Thanks again for the assistance, much appreciated.

Best regards, Spence

lammps_small_example_ALA_peptide.tar.gz (938 KB)

4

Hi Spence, thanks for checking. I’m not aware of issues between Colvars and ReaxFF, people have been using it to e.g. compute potentials of mean force of chemical reactions.

At this point I’m more inclined to look for a possible issue with forces being sent at all. Can you disable the extended mass, and just add a simple harmonic bias to the variable?

Giacomo

5

Hi Giacomo,

Thanks again for the assistance.

As requested, I disabled the extended mass and added a 2 kcal/mol/° ^2 restraint to the dihedral centered at 45°. I also included some additional simulation experiments to try and nail down the issue. In the attached input/output, it is apparent that the restraint worked and the dihedral fluctuated around 45°, while the control simulation (no harmonic bias) exhibited different behavior.

N.B., the following example input/output were a bit large, so I have sent them to Giacomo Fiorin directly. If anyone else would like them, please let me know and I can send them.

Therefore, using harmonic bias, force was sent to the atomic variables. Here is a table of additional examples I have included to help troubleshoot:

  1. Extended Mass (no) / Langevin (no) / Harmonic Bias (no): control simulation… just plain MD.

  2. Extended Mass (yes) / Langevin (yes) / Harmonic Bias (no): no forces seem to be sent to atomic dof… trajectory identical to #1.

  3. Extended Mass (no ) / Langevin (no ) / Harmonic Bias (yes): forces OK, able to restrain dihedral to target.

  4. Extended Mass (yes ) / Langevin (no ) / Harmonic Bias (yes): forces OK, able to restrain dihedral to target, fluctuations are different than without extended mass, so the extended mass appears to be OK.

  5. Extended Mass (yes ) / Langevin (yes ) / Harmonic Bias (yes): forces correct, able to restrain dihedral to target, fluctuations are different than #4, so the extended mass plus langevin appears OK.

If I had to guess (without looking at the source very closely), it may be that unless a bias is implemented, no forces are sent to the atomic dofs. I.e., with just langevin and extended mass, the function to feedback the forces from the fluctuating mass to the atomic dofs is not called.

I hope this helps… thanks again.

Best regards, Spence

6

Hi Spencer, thanks for checking with an alternative example and comparing the various options.

The result of the regular (non-extended) variable tells us that biasing forces are applied correctly to ReaxFF atoms, so that part is OK.

What I notice in your extended-Lagrangian inputs is that you keep using a damping parameter of 1 ps^{-1}, which is the default value, although provided explicitly. You may want to increase it to have a stronger temperature coupling in short tests.

Anyway, even if you run the extended mass without a thermostat, you have either left the extendedFluctuation parameter as 1 or set it to about 0.3 in either examples. That parameter, together with extendedTimeConstant, will define the mass, but again the mass may end up being too small, i.e. the extended variable will oscillate quickly around the “real” variable, and the latter won’t feel a strong influence.

I ran some tests by changing the above parameters on your system, and I do see divergent trajectories, albeit on a slow time scale.

Side note: you are not using wall potentials + extended Lagrangian for your case, but if you do please consider updating the code, because a bug was fixed recently for that particular combination (does not affect variables that do not use the extended-Lagrangian formalism).

Giacomo

7

Hi Giacomo,

Can you provide the specific values of:

extendedFluctuation

extendedTimeConstant

extendedTemperature

extendedLangevinDamping

for which you observed divergent trajectories? I would like to try and reproduce your results with my lammps / user-colvars build. I have tried (what I thought) were a wide range of these parameters, without success (i.e., the trajectories are identical to at least 1 ns of simulation time).

P.S. I have started a simulation with these values, based on your note below… it will run for 10 ns. I will start another using what you found works.

extendedFluctuation 1

extendedTimeConstant 100,000

extendedTemperature 5,032

extendedLangevinDamping 100

… those values should have a spring constant of 10 kcal/mol/°^2, and a quite large fictitious mass (with strange units since the CV has units of degree and not length).

Thanks again.

Best regards, Spence