Gradually reduce the box volume accessible to solute

[Marcelo D Poleto]

Hi,

The subject is odd, but our lab is studying aggregation patterns of 4-6 peptides in a large box of solvent. We would like to promote several aggregation events and we want to gradually reduce the volume to accessible to the peptides over time.

We are trying to design a CustomExternalForce that restrain the peptide atoms to a certain distance from the box center, as shown below. However, we are struggling with defining such a force that gradually varies.

Code: Select all

external = CustomExternalForce('K*max(0, r-dist)^2; r=sqrt(x*x+y*y+z*z)')
external.setName("Flat-external-force")
external.addGlobalParameter("K", 25 * (kilojoule_per_mole/nanometer**2))
external.addGlobalParameter("dist", 2 * angstroms)
system.addForce(external)

Does anyone have ideas or suggestions?
Thanks!

[Peter Eastman]

What you have there looks reasonable. You can gradually decrease "dist" as the simulation runs.

Code: Select all

for i in range(1000):
  simulation.context.setParameter("dist", maxDist-i*(maxDist-minDist)/1000)
  simulation.step(100)

Is there something specific that isn't working?

[Marcelo D Poleto]

Hi Peter,

Your suggestion makes total sense. I will try to implement it and I will come back with the final structure. Could you explain a bit more the rationale behind the expression "maxDist-i*(maxDist-minDist)/1000"?

Also, we will be working with a periodic system not centered at (0,0,0) and the implementation above was allowing the peptides to cross the PBC, which messes up with the protocol principle. So I tried to center the potential at the boxcenter by offsetting the coordinates by a boxcenter value, but the peptides are just constantly being "diagonally pulled" across the box.

Code: Select all

external = CustomExternalForce("K*max(0, r-dist)^2; r=sqrt(a*a+b*b+c*c); \
                               a=(x-centerx); b=(y-centery); c=(z-centerz)")
external.setName(flat_name)

external.addGlobalParameter("centerx", boxcenters[0].value_in_unit(angstrom))
external.addGlobalParameter("centery", boxcenters[1].value_in_unit(angstrom))
external.addGlobalParameter("centerz", boxcenters[2].value_in_unit(angstrom))
external.addGlobalParameter("K", flat_constant)
external.addGlobalParameter("dist", 10 * angstroms)

Do you have an idea of what I might be missing here?

[Peter Eastman]

That expression just linearly reduces the distance. When i is 0 it equals maxDist. When i is 1000 it equals minDist. You can use whatever expression you want for your protocol.

To compute the displacement while using periodic boundary conditions, compute it as r=periodicdistance(x, y, z, centerx, centery, centerz).

[Marcelo D Poleto]

Hi Peter,

Thanks for the clarification. It worked very well! I used pretty much all of your suggestions here, but for future reference this is how my code looks like:

Code: Select all

flat_vol_constant = 100 * (kilojoule_per_mole/nanometer**2)
flat_name = "Flat-external-force"

external = CustomExternalForce("Kvol*max(0, r-dist)^2 ; \
                        r=periodicdistance(x, y, z, centerx, centery, centerz)")

external.setName(flat_name)
external.addGlobalParameter("centerx", 0.0*angstrom) #aggregation at origin
external.addGlobalParameter("centery", 0.0*angstrom) #aggregation at origin
external.addGlobalParameter("centerz", 0.0*angstrom) #aggregation at origin
external.addGlobalParameter("Kvol", flat_vol_constant)
external.addGlobalParameter("dist", maxdist) # start allowing maximum distance
system.addForce(external)

u_tmp = MDAnalysis.Universe(pdbfile)
vol_selection = u_tmp.select_atoms("resname POT")
for i in vol_selection.indices:
    external.addParticle(i, [])

And later, after defining the simulation object:

Code: Select all

for s in range(decrease_stages):
        accessible_dist = maxdist - s*(maxdist-mindist)/decrease_stages
        simulation.context.setParameter("dist", accessible_dist)
        simulation.step(nsteps_decrease_per_stage)

Thank you for the insights!