Page 1 of 1
Gradually reduce the box volume accessible to solute
Posted: Fri Nov 03, 2023 8:45 am
by mdpoleto
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!
Re: Gradually reduce the box volume accessible to solute
Posted: Fri Nov 03, 2023 10:32 am
by peastman
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?
Re: Gradually reduce the box volume accessible to solute
Posted: Sat Nov 04, 2023 10:04 am
by mdpoleto
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?
Re: Gradually reduce the box volume accessible to solute
Posted: Sat Nov 04, 2023 11:19 am
by peastman
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).
Re: Gradually reduce the box volume accessible to solute
Posted: Tue Nov 07, 2023 1:12 pm
by mdpoleto
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!