The issue is that when running NPT w/ restraints, the restraint coordinates are scaled along with the regular coordinates (by molecule by default). Also, by default, LEaP adds too little water so upon pressure scaling, the box has to shrink and therefore the restraint coordinates shrink. No matter your force constant, you cannot get back to the original since the refc was scaled too... You could get around this by doing short MD runs where the -refc at each restrt is always the original coords (and not the restrt). This never bit us in the old days since our equilibration protocols were not as elaborate (or sensible) as they are now.
The hack I used in the past (which was documented in the very old DNA tutorial) was to edit the prmtop to make the dimer a single molecule. This is a bit tricky to edit by hand as multiple changes are required in the prmtop. This could be added to parmed...
Another alternative is during the solvateOct or solvateBox to set the closeness to 0.97; this leads to a more realistic density. refc coordinates may still be scaled, but only modestly.
The final alternative is to comment out the scaling of the refc coordinates (in the sander and pmemd codes). I've done this and seen no ill effects, but perhaps my group should validate and fix this finally rather than letting the issue languish. I do not know the origin of the decision to scale the refc coordinates in addition to the coords, but it's been there from the beginning.
A number of the AMBER developers have our (final day of virtual) annual meeting tomorrow and I'll bring it up. To date I've hacked around this rather than solving the problem, so for that I am sorry this bit you and other users.
--tec3
\-/ Professor Thomas E. Cheatham, III
-/- Department of Medicinal Chemistry, College of Pharmacy
/-\ Director, Research Computing and CHPC, UIT, U of Utah
\-/
-/- tec3.utah.edu
http://www.chpc.utah.edu/~cheatham
/-\ SKH-4914 (801) 587-9652; FAX: (801) 585-6208
\-/ INSCC-410 (801) 585-6318; FAX: (801) 585-5366
________________________________________
From: David A Case <david.case.rutgers.edu>
Sent: Thursday, February 25, 2021 7:02:24 PM
To: AMBER Mailing List
Subject: Re: [AMBER] dimer shifts under NPT ensemble
On Fri, Feb 26, 2021, Qinghua Liao wrote:
>
>My understanding is that when the ensemble switched from NVT (heating)
>to NPT (equilibration), the volume of the water box shrinked,
>which pushed the two monomers getting closer. I checked the box
>dimensions, each dimension shrinked by 2 angstrom.
>
>Even though the shift is just a little, I am still worried that this
>shift could still somehow mess up the interactions at the interface,
>leading to a relatively unstable dimer.
I don't know of any easy way to avoid this problem altogether. Don't use
really strong restraints in NPT steps, to allow the dimer interface to relax
as the box shrinks.
>
>I also tried to reduce the closeness of waters (0.75 angstrom) when
>generating the topology and coordinate files, as it could add more water,
>and the water box might shrink less when switching from NVT to NPT
>ensembles. But it was still the same.
When you say "it was still the same", do you mean that nothing at all
changed? Using a closeness of 0.75 should increase the initial density,
and hence reduce the amount by which the box needs to shrink. If nothing
changes when you change the closeness parameter, maybe there is something
wrong at the tleap step.
Finally, do you know that the dimer is unstable under your simulations, or
is that just a "worry"?
....dac
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Feb 25 2021 - 20:00:01 PST
