Also, isn't cut_factor (pbc.F90) supposed to account for this? Looks like
the ucell matrix remains upper diagonal under all circumstance?
114 if (is_orthog .ne. 0) then
115 ucell(:, :) = 0.d0
116 ucell(1, 1) = a
117 ucell(2, 2) = b
118 ucell(3, 3) = c
119 cut_factor(:) = 1.d0
120 else
121 factor = PI / 180.d0
122 ucell(1, 1) = a
123 ucell(2, 1) = 0.d0
124 ucell(3, 1) = 0.d0
125 ucell(1, 2) = b * cos(factor * gamma)
126 ucell(2, 2) = b * sin(factor * gamma)
127 ucell(3, 2) = 0.d0
128 ucell(1, 3) = c * cos(factor * beta)
129 ucell(2, 3) = (b * c * cos(factor * alpha) - ucell(1, 3) * &
130 ucell(1, 2))/ucell(2, 2)
131 ucell(3, 3) = sqrt(c * c - ucell(1, 3) * ucell(1, 3) - ucell(2,
3) * &
132 ucell(2, 3))
133
134 ! Cut factors are used to correct for "spherical cutoff
protrusion" into
135 ! adjacent unit cells. The problem is that the point at which
a cutoff
136 ! sphere is tangent to a unit cell side is not the contact
point for
137 ! projection of an orthogonal x, y, or z vector in a
nonorthogonal unit
138 ! cell. We thus have to increase the cutoff a bit to allow for
the longer
139 ! distance for the orthogonal projection.
140
141 cut_factor(1) = 1.d0 / (sin(factor * beta) * sin(factor *
gamma))
142 cut_factor(2) = 1.d0 / (sin(factor * alpha) * sin(factor *
gamma))
143 cut_factor(3) = 1.d0 / (sin(factor * alpha) * sin(factor *
beta))
144 end if
On Fri, Mar 18, 2022 at 12:21 PM Scott Le Grand <varelse2005.gmail.com>
wrote:
> Finally, the fix is surprisingly atomic modulo working the fixed point
> math into cellOfffset.
>
> On Fri, Mar 18, 2022 at 12:20 PM Scott Le Grand <varelse2005.gmail.com>
> wrote:
>
>> Also the fix is in(tm), but I need to address gpu_update_box and I need
>> to know if ucell matrices are still upper diagonal because if not that will
>> require some further work as I had to update the cell offsets to be
>> precalculated to make fixed point consistent and properly bounded so there
>> can't ever be an overflow.
>>
>>
>>
>> On Fri, Mar 18, 2022 at 12:17 PM Scott Le Grand <varelse2005.gmail.com>
>> wrote:
>>
>>> First up, thanks for replying! Appreciated!
>>>
>>> But... For triclinic cells, I never saw a ucell matrix full occupied.
>>> u(2,1), u(3,1), and u(3,2) are always zero, yes?
>>>
>>> If so, then the column space defines a, b, and c by the respective
>>> column lengths that *are* the distances between those separating planes.
>>> Check it out for yourself. Can you point me to a case that violates this?
>>> If so, we have to fix the whole shebang, not just this routine. Because I
>>> agree with you if we have fully occupied ucell matrices, otherwise no.
>>>
>>>
>>>
>>> On Fri, Mar 18, 2022 at 12:10 PM David Cerutti <dscerutti.gmail.com>
>>> wrote:
>>>
>>>> What's going on here is that not all unit cells are orthorhombic--such
>>>> would imply that the box angles alpha, beta, and gamma are always 90
>>>> degrees. The general unit cell is triclinic, where alpha, beta, and
>>>> gamma
>>>> can take on any values (with some degree of excluded territory, as
>>>> certain
>>>> combinations would fail to make a shape that can exist). The code you
>>>> are
>>>> seeing are seeing was created in response to another aspect of the
>>>> original
>>>> pmemd.cuda pair list that, if I recall correctly, was missing
>>>> interactions
>>>> from time to time due to a miscommunication between the way the pair
>>>> list
>>>> was generated versus the way the box had changed. I don't recall if it
>>>> was
>>>> creating energy drift or causing an outright crash. It was also around
>>>> the
>>>> time I added the -AllowSmallBox flag: note that the code will not allow
>>>> users to run simulations that are less than three pair list subdivisions
>>>> along any dimension without an explicit instruction to do so, as I
>>>> wasn't
>>>> sure that the way the particles get re-imaged at the beginning of each
>>>> time
>>>> was safe if there are only two subdivisions.
>>>>
>>>> But the HessianNorms function is calculating the actual thickness of
>>>> each
>>>> subdomain cell, as defined by the distance separating the planes of
>>>> opposing faces. This is integral to the reason that, unless you do some
>>>> sophisticated stuff like the honeycomb idea I had, a truncated
>>>> octahedron
>>>> is not as efficient as the 22% reduction in overall particle count would
>>>> imply: the total volume of each cell needs to expand as it tilts, and
>>>> this
>>>> implies more tiles to cover the same volume. It's still advantageous to
>>>> use a truncated octahedron, as the Hilbert Space curve does pretty well
>>>> to
>>>> localize what is needed and of course the cutoff shell of any given
>>>> particle remains spherical.
>>>>
>>>> What's the ETA on the energy drift fix? Is this an atomic change, or
>>>> does
>>>> it also come with other refactoring / rollback of various precision
>>>> trims?
>>>>
>>>> Dave
>>>>
>>>>
>>>> On Fri, Mar 18, 2022 at 2:14 PM Scott Le Grand <varelse2005.gmail.com>
>>>> wrote:
>>>>
>>>> > Let's look at this, here's the new code:
>>>> > HessianNorms(invU, cdepth);
>>>> > double skin = cdepth[0]/gpu->sim.xcells - gpu->sim.cut;
>>>> > double yskin = cdepth[1]/gpu->sim.ycells - gpu->sim.cut;
>>>> > double zskin = cdepth[2]/gpu->sim.zcells - gpu->sim.cut;
>>>> >
>>>> > The array cdepth is filed with the new a, b, and c by means of an
>>>> external
>>>> > call. But as far as I know, all unit cells are orthorhombic, so we've
>>>> > always been able to use:
>>>> > // Recalculate a, b, c
>>>> > double a =
>>>> > pNTPData->ucell[0];
>>>> > double b =
>>>> > sqrt(pNTPData->ucell[1] * pNTPData->ucell[1] +
>>>> >
>>>> > pNTPData->ucell[4] * pNTPData->ucell[4]);
>>>> > double c =
>>>> > sqrt(pNTPData->ucell[2] * pNTPData->ucell[2] +
>>>> >
>>>> > pNTPData->ucell[5] * pNTPData->ucell[5] +
>>>> >
>>>> > pNTPData->ucell[8] * pNTPData->ucell[8]);
>>>> >
>>>> > because of the constraints that imposes on the column space of the
>>>> unit
>>>> > cell matrix.
>>>> >
>>>> > So either the first code clause is overkill and repeated work or we
>>>> have to
>>>> > revisit the previous orthorhombic assumption. Which is it?
>>>> >
>>>> > But even then, why would one change the unit cell matrix without
>>>> updating
>>>> > the a, b, and c derived from it? And yes, I'm impatient w/r to fixing
>>>> bugs.
>>>> > Guilty.
>>>> >
>>>> >
>>>> > On Fri, Mar 18, 2022 at 10:42 AM Scott Le Grand <
>>>> varelse2005.gmail.com>
>>>> > wrote:
>>>> >
>>>> > > If the unit cell changes, wouldn't a, b, and c change as well since
>>>> they
>>>> > > are derived from said unit cell?
>>>> > >
>>>> > > As for the speed of response. If it's important to anyone, they'll
>>>> > > respond. If not, well they can make their principled stand against
>>>> me and
>>>> > > my impatience later. Otherwise, I acknowledge your point but also my
>>>> > > okeydokey here, nobody is debugging this thing. It's an easy fix.
>>>> And I
>>>> > > need to do that fix to make it exploit the fixed point local
>>>> coordinate
>>>> > > system, but I could leave the seemingly sloppy incorrect behavior
>>>> here.
>>>> > > Driveby contributions to PMEMD like this however are just not good
>>>> IMO
>>>> > and
>>>> > > why I threw my hands up in frustration after trying to refactor for
>>>> > > performance in 2020. It's a Sisyphusian effort and the tsunami of
>>>> > > unattended warning messages and dead code in the thing alone is
>>>> enough to
>>>> > > make one run screaming into the night IMO. It was a single
>>>> unitialized
>>>> > > variable that broke PMEMD for Silicon Therapeutics in 2020. But
>>>> also, why
>>>> > > is there is so much dead code?
>>>> > >
>>>> > > The author of this code added a new and completely superfluous
>>>> method of
>>>> > > recalculating a, b,and c without exploiting the properties of the
>>>> > > orthorhombic coordinate system to do it more simply as it is done in
>>>> > > gpu_pressure_scale, which itself does the same thing internally as
>>>> > > gpu_update_box if there's an MC barostat. And I hope you can
>>>> forgive me
>>>> > for
>>>> > > being annoyed at the needless repetition of functionality here w/o
>>>> any
>>>> > > notes as to why this was necessary.
>>>> > >
>>>> > > On Fri, Mar 18, 2022 at 10:22 AM David A Case <
>>>> david.case.rutgers.edu>
>>>> > > wrote:
>>>> > >
>>>> > >> On Fri, Mar 18, 2022, Scott Le Grand wrote:
>>>> > >>
>>>> > >> >Okeydokey, nobody cares. Therefore it's a bug and I'm going to
>>>> fix it.
>>>> > >> Test
>>>> > >> >will be updated to reflect this. What SM revision should I use
>>>> for the
>>>> > >> test
>>>> > >> >reference files?
>>>> > >>
>>>> > >> ??!? How fast do you expect people to respond? I looked at this
>>>> and
>>>> > >> didn't
>>>> > >> see any obvious error -- the two routines are called at different
>>>> > places:
>>>> > >> gpu_pressure_scale only when ntp>0 and barostat=1, whereas
>>>> > gpu_update_box
>>>> > >> is only called when re-weighting trajectories.
>>>> > >>
>>>> > >> But I don't know the code here. Clearly, you must be convinced
>>>> that
>>>> > a,b,c
>>>> > >> need to be updated when reweight is > 0. I'll try to look at what
>>>> > >> behavior
>>>> > >> changes when you make the change. But I'm travelling to the ACS
>>>> meeting
>>>> > >> tomorrow, so you may disappointed if you require immediate
>>>> feedback.
>>>> > >>
>>>> > >> ....dac
>>>> > >>
>>>> > >>
>>>> > >> _______________________________________________
>>>> > >> AMBER-Developers mailing list
>>>> > >> AMBER-Developers.ambermd.org
>>>> > >> http://lists.ambermd.org/mailman/listinfo/amber-developers
>>>> > >>
>>>> > >
>>>> > _______________________________________________
>>>> > AMBER-Developers mailing list
>>>> > AMBER-Developers.ambermd.org
>>>> > http://lists.ambermd.org/mailman/listinfo/amber-developers
>>>> >
>>>> _______________________________________________
>>>> AMBER-Developers mailing list
>>>> AMBER-Developers.ambermd.org
>>>> http://lists.ambermd.org/mailman/listinfo/amber-developers
>>>>
>>>
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Fri Mar 18 2022 - 13:00:02 PDT
