Hi Tom,
that's not quite the issue I'm talking about, which definitely goes away
when you use iwrap=0.
I think the pressure scaling issue only happens with multiple solute
molecules, right? I don't think we scale atom positions (by default
anyway), just shift entire molecules. So I have only seen the problem you
mention for the helix when you have a multimer or complex, and they shift
relative to each other (and violate restraints on restart). That one is
also bad but doesn't usually cause thousands of kcal restraint energy. A
possible easy fix for the constant P scaling with restraints would be to
have parmed or something make sure that the restrained atoms are all in the
same "molecule" in the prmtop, or change the pressure scaling code to only
change things past NSPSOL, for example... we edit this when needed to avoid
the pressure issues. I also agree with not scaling refc, that's
probably the better physics since your reference should never change, IMO.
On Mon, Aug 8, 2022 at 11:00 PM Thomas Cheatham <tec3.utah.edu> wrote:
>
> > I've run into a problem (not technically a bug, but a set of input
> > combinations that leads to inability to restart a run). I'm emailing to
> > discuss how best to fix it, or maybe someone knows what I'm doing wrong.
>
> It was a choice, and one well before my time, that I do not understand...
>
> It does not have to do with iwrap, netcdf, ASCII, or any of that, I think.
> Less of a problem with molecule pressure scaling than atomic positional
> scaling, depending on restraints (since normally you do not apply
> restraints across multiple molecules).
>
> But, a choice was made a long time ago to scale for pressure both the
> coordinates and the refc coordinates, so your refc changes if your box
> (with defaults) is not at the correct density. Two fixes, get the box
> closer to correct density at start with the setbox I think closer or
> something else to 0.98 or comment out the code that scales the refc
> coordinate scaling. The box scaling may shift the relatives of the
> restraint coordinates, but they should come back.
>
> With current implementation, say you were restraining an alpha helix and
> box had to shrink a fair bit, the helix would shrink to match the
> restraint coordinates since they were shrinking and then upon restart they
> would suddenly be well away from the original unscaled refc.
>
> I've talked about this on forum (10+ years back) but did not fix (my bad)
> or put in option to not scale the refc coordinates (which I feel is
> correct). On a step by step basis, if refc is constant it will put the
> restrained atoms back w/out huge forces or issues.
>
> --tec3
>
> (yes, may be worse with iwrap, but issue is present regardless unless I or
> someone else did fix and I've forgotten about it)
>
> > my problem is that when you set iwrap=1, coordinates can get changed even
> > for atoms with positional restraints. The coordinates of these atoms can
> be
> > translated substantially in the restart compared to the input
> coordinates.
> > Then, when you try to restart, the new coordinates no longer match the
> > reference coordinates (refc file), and the simulation blows up with
> > huge restraint energy. There is no simple way to fix it and successfully
> > continue the previous run. It's not a problem during the MD, only when
> you
> > try to restart. For me, using iwrap=1 quickly fails during a multi-step
> > equilibration (sorry, "relaxation"). iwrap=0 is fine, but not good for
> very
> > long MD.
> >
> > - you can use the restart file of the previous run as the inpcrd AND the
> > refc. This works fine now, but I've never really liked the idea of
> > continually changing what you call your "reference" structrure.
> > - we could enforce iwrap=0 when using ntr=1 (but this can lead to
> > coordinate overflow, especially for solvent)
> > - we could shift the reference coordinates and save them to a "reference
> > restart" file? seems too complicated.
> > - we could skip translating atoms that have restraints applied? they
> aren't
> > moving far anyway.
> >
> > other ideas?
> >
> > the relevant code (I think) is in pbd.F90, wrap_molecules()
> >
> > do i = offset + 1, offset + mol_atm_cnt
> > atm_id = gbl_mol_atms_lists(i)
> > crd(1, atm_id) = crd(1, atm_id) + tran(1)
> > crd(2, atm_id) = crd(2, atm_id) + tran(2)
> > crd(3, atm_id) = crd(3, atm_id) + tran(3)
> > end do
> >
> >
> >
> > other ideas?
> > thanks
> > carlos
> > _______________________________________________
> > AMBER-Developers mailing list
> > AMBER-Developers.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber-developers
> >
>
> \-/ 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
>
>
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Tue Aug 09 2022 - 06:30:03 PDT
