Atom Stuck at Harmonic Wall

Hi all,

I am simulating the dissolution of Calcium off a mineral, with respect to rare event observations. Because of that, I would like the atom to detach and reattach without drifting off into the water and never coming back. To do such, I have implemented 1 dimensional Metadynamics and ABF. From this, I have used the Harmonic Walls command in Colvars and as a Lammps region (separately) with varying force constants and parameters. I have also put force constant restraints on the custom function for my defined Colvar as well. The goal is to have a spherical region around it.

However, after running the simulation for a few nS the atom interacts with the “wall” and appears to get stuck at that position. The pressure increases, all atoms in the simulation lock in position and barley vibrate, even though the simulation still runs.

I am curious if anyone has some idea or reason to why the atom interacts with the wall, then gets stuck and “freezes” at its position not being able to break off the wall.

The input files have the extra commented out sections reflecting the different methods used, but not with all the different values and trials.

Lastly, the goal is not to define a Colvar based on a distance, rather I am trying to look at and use the coordination number.

Thank you all in advance,
Tyler

LAMMPS INPUT:

CaCO3

units metal

boundary p p f

atom_style full
bond_style harmonic
angle_style class2
improper_style distance
pair_style hybrid/overlay coul/long 9. lj/cut 9. &
buck/mdf 6.0 9.0 lj/mdf 6.0 9.0

read_data restartrelax1.data
#read_restart lammps.restart.5000

reset_timestep 0

group slab type 1 2
group OC type 3
group OW type 4
group CA id 254
group2ndx groups.ndx OC OW CA

BONDS

bond_coeff 1 20.424650 1.3042000
bond_coeff 2 22.965000 1.0120000

ANGLES

angle_coeff 1 120.000 6.617 0.0 0.0
angle_coeff 1 bb 12.818 1.3042 1.3042
angle_coeff 1 ba 1.53319 1.53319 1.3042 1.3042

angle_coeff 2 113.240 1.64568 0.0 0.0
angle_coeff 2 bb 0.0 0.0 0.0
angle_coeff 2 ba 0.0 0.0 0.0 0.0

IMPROPERS

improper_coeff 1 13.647 360.00

PAIR COEFFS

pair_coeff * * coul/long
pair_coeff 1 3 buck/mdf 3161.6335 0.271511 0.0 6.0 9.0
pair_coeff 1 4 lj/mdf 0.00095 3.35 6.0 9.0
pair_coeff 3 3 buck/mdf 63840.199 0.198913 27.89901 6.0 9.0
pair_coeff 3 4 buck/mdf 12534.455133 0.202 0.0 6.0 9.0
pair_coeff 3 5 buck/mdf 340.0 0.217 0.0 6.0 9.0
pair_coeff 4 4 lj/cut 0.0067400000 3.165492
comm_modify vel yes
pair_modify tail yes

kspace_style pppm 1.0e-5
kspace_modify slab 3.0
thermo 100
thermo_style multi

dump 1 all xtc 100 lammps.xtc
#dump_modify 1 unwrap yes

#creating shpere region
#region mySphere sphere -2.958009 0.346599 -10.626204 7.75 side in

name type x y z rad inside

timestep 0.001
fix 1 all nvt temp 298.15 298.15 0.1 tchain 5
fix 2 slab recenter INIT INIT INIT shift all units box
fix 3 all wall/reflect zlo -70.0 zhi 70.0
#new wall fix
#fix 4 CA wall/region mySphere harmonic 1000.0 0.0 4.5

sig eps cut

fix 4 all colvars COLVAR_FILE #input rest.colvars.state

#minimize 1e-4 1e-4 1000 1000
run 6000000

write_data restart2.data nocoeff

COLVARS INPUT:

collective variable example: monitor distances

colvarsTrajFrequency 500 # output every 500 steps #(0.1 psec)
colvarsRestartFrequency 10000

indexFile groups.ndx

colvar {
name one
customFunction (y + 1) / ( x )
lowerBoundary 0.0
upperBoundary 10.0
width 0.05
extendedLagrangian yes
extendedTemp 298
extendedFluctuation 0.05
extendedTimeConstant 0.04
coordNum {
name y
cutoff 3.3
expNumer 10.0
expDenom 64.0
group1 {indexGroup CA }
group2 {indexGroup OC }
}
coordNum {
name x
cutoff 4.0
expNumer 8.0
expDenom 40.0
group1 {indexGroup CA }
group2 {indexGroup OW }
}
}

colvar {
name wall
lowerBoundary 0.0
upperBoundary 10.0
distance {
group1 {indexGroup CA }
group2 {dummyAtom (-2.958009, 0.346599, -10.626204) }
}
}

harmonicWalls {
name harm_wall
colvars wall
upperWalls 5.0
upperWallConstant 1.0
}

abf {
name abf_dist
colvars one
fullSamples 100
outputFreq 500
historyFreq 500
}

#metadynamics {

name meta_dist

colvars one

hillWeight 0.025

writeFreeEnergyFile yes

keepFreeEnergyFiles on

newHillFrequency 2500

hillWidth 2.5

keepHills

wellTempered on

biasTemperature 2980

#}

This is likely a Q for the Colvars folks (CCd).

Steve

Hi Tyler, sorry I missed this (and thanks, Steve!).

There are a few things that may be misbehaving. The first to look at is how the variable depends on the atomic coordinates: you are using a non-linear function (y+1)/x of two non-linear functions (coordNum), so it is difficult for me to get an idea of the typical range of values.

One main feature of a coordNum is that its atomic gradients become vanishingly small when you approach either zero or the maximum of the function. This will affect how easily the atoms respond to forces applied to the collective variable. Since you are using a harmonic wall anyway, can you elaborate on why a distance-based variable wouldn’t work for you?

There are some subtle issues of metadynamics and walls or boundaries, summarized here:
https://colvars.github.io/colvars-refman-lammps/colvars-refman-lammps.html#sec:colvarbias_meta_boundaries
(They may not apply to your case, though, because your wall is within your grid boundary.)

ABF, but more specifically extended-system ABF may suffer from some inaccuracies near walls, and there is the option of applying the wall force directly to the atoms instead of the virtual mass (new keyword “bypassExtendedLagrangian”):
https://github.com/Colvars/colvars/pull/299

For the latter feature especially you will probably need to take the latest LAMMPS snapshot and update it with the latest Colvars snapshot. This is also more conveniently downloadable as a Git branch here:

https://github.com/giacomofiorin/lammps/tree/colvars-update-2020-03-17
(the branch is automatically generated each day and then removed; update the date to today’s date if needed).

Hope it helps. Stay safe and healthy.

Giacomo

Hey all,

Thank you for the responses! Always helpful! I hope everyone is doing well with whole Corona situation…

(Apologies for my linear and listed response… Can be hard responding to multiple threads at once, and I hope it doesn’t come across as rude)

Giacomo…

With respect to…

  1. non-linear function: This function started out as a ratio of the two coordination numbers. However when the reaction progresses “y”: coordination to carbonate oxygen’s, eventually will approach and could be zero. Remembering correctly, this would cause funky things in the calculations, so a +1 was added… However this may not address what you were referencing… ie how the atomic gradients are working. Anyways a range of values: ~0.2 for a detached calcium ion, and ~1.2-1.6 attached at the site. Most common values are the ones between the two previously mentioned.

  2. Why distance won’t cut it…: The goal of the project to initially get a sense of the energy landscape in 1 and 2 dimensions, and use this to guide our Forward Flux Sampling calculations to derive dissolution rates based on coordination numbers (not as a summation of “Boolean” (kinda) 1 or 0, coordinating or not) and not distance. As proposed, there are several advantages to this as compared to distance. To briefly state. Furthering that, you stated “This will affect how easily the atoms respond to forces applied to the collective variable.” Could you go in more depth on this? Or at maybe point me in the direction of where to read about this? So I can fully understand what could be going on?

  3. I will check out those Metadynamics links you shared…

  4. I have messed around with ABF. Where I haven’t had the best success in it, when generating “acceptable” energy landscapes, I have not seen the bypassExtendedLagrangian function yet. Thank you for that input and I will check it out!

  5. I will check out the snapshot as well!

Furthering this topic… In the last month or so I came across the Drag fix. This shed some light on what was going on. I believe that as the energy landscape progresses, the values near the extremes of the colvar increase. Specifically the detached state. I have noticed that where the wall may work for the first interaction, then it doesn’t work for following ones. I am thinking that at the first interaction that the force in the energy landscape is lower than the force constant on the wall. So it works. Then the force from the energy landscape gets greater than the force from the wall, “trapping” the ion in between the wall and the force from the energy landscape. Where you can’t change the force constant between runs in Lammps, you can change the fixes. So, when I noticed that the atom got stuck, I would repeat that specific run with a higher force constant on the drag fix. Hopefully, I can hone in on the proper force constant for the entirety being able to use the harmonic wall in the Colvars system because I think that is more representative of real world forces than others.

Hope all is well and everything made sense,

Tyler Schmidt

Hi Tyler, a highly non-linear function xi® with sigmoidal-like shape (like that used for coordNum) has derivative dxi®/dr that becomes very small for r very small or r very large. The magnitude of dxi®/dr is largest near the inflection point of the curve. Any force applied to xi® will be multiplied (due to the chain rule) by dxi®/dr before being transformed into an atomic force, thus it will vary greatly in magnitude.

When you use a coordNum variable you have a sum over many such functions, and at any given time frame not all of them may be near their inflection point. This means that there will be several occasions where none of these terms will be at the inflection point, and the variable does not increase/decrease easily. This is generally not a problem if you are applying only small forces to the variable (e.g. to compute a PMF), because although the effect of the applied force is negligible, diffusion will eventually move the system to a new configuration. But applying a wall force to fix a certain configuration may not have the result that’s intended.

About the use of fix drag: without knowing the parameters you chose, keep in mind that it’s best not to apply non-conservative forces within the range of the PMF. At least, do not consider those points in the PMF where fix drag is acting.

Giacomo