Energy minimization with fix rigid/small

[demkowicz]

Under the "Restrictions" heading on the bottom of the minimize command manual entry (https://docs.lammps.org/minimize.html), there is a statement you'd be interested in implementing energy minimization under the fix rigid constraint. I believe this could be achieved pretty easily, at least for fix rigid/small using derivative methods such as conjugate gradient. I'm not familiar with LAMMPS source code, but, if there is interest, I would be happy to work with someone who knows it to implement this functionality. I would personally very much like to see such functionality in LAMMPS.
M. J. Demkowicz (demkowicz@tamu.edu)

[athomps]

I could certainly give you some guidance. I suggest starting out by looking at two functions: MinLineSearch::alpha_step() and Min::energy_force(). Let me know how you feel after that :-)

[sjplimp]

You should look at this paper for HOOMD. I don't have a copy at the moment, but it has a section on energy min with rigid body constraints. My recollection is it uses FIRE (damped dynamics) not a CG based method.
@ndtrung81 Trung can comment as he is the lead author. I think this is likely a better match to how LAMMPS would use MPI (or GPU) parallelism.

Nguyen, T. D., Phillips, C. L., Anderson, J. A., & Glotzer, S. C. (2011). Rigid body constraints realized in massively-parallel molecular dynamics on graphics processing units. Computer Physics Communications, 182(11), 2307-2313. https://doi.org/10.1016/j.cpc.2011.06.005

[ndtrung81]

Correct, in that paper we extended the FIRE algorithm (E. Bitzek et al, Phys. Rev. Lett. 97 (2006)) to account for the contribution of the body rotational degrees of freedom to the power (body torque times angular velocity in the space-fixed frame). I think it would be possible to either extend the current MinFire to handle rigid bodies (when fix rigid is present), or to have a separate MinFireRigid deriving from MinFire. The former approach would allow energy minimization of systems composed of non-rigid parts and rigid bodies but need changes to MinFire. In the latter, it is difficult to ensure two MinFire* instances update the atoms in the non-rigid parts and those in the bodies consistently.

[demkowicz]

Aidan, I had a look, but it doesn’t really tell my anything. I’ve written my own implementations of CG minimization, so I understand the algorithm. However, I’ve never worked with LAMMPS source code. — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.******@***.***> / 979-458-9845 / RDMC 212 On May 31, 2023, at 12:20 PM, Aidan Thompson ***@***.******@***.***>> wrote: I could certainly give you some guidance. I suggest starting out by looking at two functions: MinLineSearch: : alpha_step() and Min: : energy_force(). Let me know how you feel after that :-) — Reply to this email directly, view it on GitHub, ZjQcmQRYFpfptBannerStart This Message Is From an External Sender This message came from outside your organization. ZjQcmQRYFpfptBannerEnd I could certainly give you some guidance. I suggest starting out by looking at two functions: MinLineSearch::alpha_step() and Min::energy_force(). Let me know how you feel after that :-) — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/lammps/lammps/issues/3801*issuecomment-1570617700__;Iw!!KwNVnqRv!CS1-3qr-WUSpGAHhvS46Z0_M2ggi3JjpVOZrUwENB62OvL7Rd174GRLcqSvMoOizQLnQYu7gT1ZS4opPyEsmcmyeJ4A$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AD3OZLXADXOWTIG577EITJTXI54WXANCNFSM6AAAAAAYVJ665Y__;!!KwNVnqRv!CS1-3qr-WUSpGAHhvS46Z0_M2ggi3JjpVOZrUwENB62OvL7Rd174GRLcqSvMoOizQLnQYu7gT1ZS4opPyEsm9H-cLGA$>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[demkowicz]

I’m not familiar with the FIRE algorithm. However, overdamped dynamics is very close to constant increment steepest descent. Glancing over the paper, I see it has standard rigid body dynamics with the only apparent novelty being in the integration scheme. This makes me think it may be possible to do a quick kluge within lammps script where you run 1 tilmestep at a time in a loop. On each step, you set the COM velocity to the COM force and the angular velocity to the torque for each molecule. Then you adjust the tilmestep to limit the max displacement and rotation. This may be a relatively simple thing to try first with no need to alter the source code. — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.******@***.***> / 979-458-9845 / RDMC 212 On May 31, 2023, at 4:07 PM, Steve Plimpton ***@***.******@***.***>> wrote: This Message Is From an External Sender This message came from outside your organization. You should look at this paper for HOOMD. I don't have a copy at the moment, but it has a section on energy min with rigid body constraints. My recollection is it uses FIRE (damped dynamics) not a CG based method. @ndtrung81<https://urldefense.com/v3/__https://github.com/ndtrung81__;!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfcY-bfDCg$> Trung can comment as he is the lead author. I think this is likely a better match to how LAMMPS would use MPI (or GPU) parallelism. Nguyen, T. D., Phillips, C. L., Anderson, J. A., & Glotzer, S. C. (2011). Rigid body constraints realized in massively-parallel molecular dynamics on graphics processing units. Computer Physics Communications, 182(11), 2307-2313. https://doi.org/10.1016/j.cpc.2011.06.005<https://urldefense.com/v3/__https://doi.org/10.1016/j.cpc.2011.06.005__;!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfcpxsTSGI$> — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/lammps/lammps/issues/3801*issuecomment-1570956487__;Iw!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfcZqjwuos$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AD3OZLUS7GDZVE3XVAKK46DXI6XH7ANCNFSM6AAAAAAYVJ665Y__;!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfc3Dy3nWE$>. You are receiving this because you authored the thread.

[demkowicz]

Trung, perhaps you would be interested in implementing this in LAMMPS? As far as I know, LAMMPS currently does not have any capability of doing energy minimization with rigid body constraints. — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.******@***.***> / 979-458-9845 / RDMC 212 On Jun 1, 2023, at 10:47 AM, Trung Nguyen ***@***.******@***.***>> wrote: This Message Is From an External Sender This message came from outside your organization. Correct, in that paper we extended the FIRE algorithm (E. Bitzek et al, Phys. Rev. Lett. 97 (2006)) to account for the contribution of the body rotational degrees of freedom to the power (body torque times angular velocity in the space-fixed frame). I think it would be possible to either extend the current MinFire to handle rigid bodies (when fix rigid is present), or to have a separate MinFireRigid deriving from MinFire. The former approach would allow energy minimization of systems composed of non-rigid parts and rigid bodies but need changes to MinFire. In the latter, it is difficult to ensure two MinFire* instances update the atoms in the non-rigid parts and those in the bodies consistently. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/lammps/lammps/issues/3801*issuecomment-1572302613__;Iw!!KwNVnqRv!C4QsJJ9ujszb8Riu4EY3-Bfzkkc0ab-TGLxU9-eWDgckw77WIh5X0dVPJhliLf-a1BiiUu_NYJoZqRc-u1aJzPIQFEI$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AD3OZLTCPLCE5AUERCBMHA3XJC2QNANCNFSM6AAAAAAYVJ665Y__;!!KwNVnqRv!C4QsJJ9ujszb8Riu4EY3-Bfzkkc0ab-TGLxU9-eWDgckw77WIh5X0dVPJhliLf-a1BiiUu_NYJoZqRc-u1aJqGLPIvE$>. You are receiving this because you authored the thread.

[demkowicz]

BTW, Steve, why do you anticipate parallelization issues for conjugate gradient? Is it because the line search might lead to large displacements that could require remaking the linked list? This is possible in my experience, but not likely. The performance advantage of CG over steepest descent is rather, uh… steep! ;) — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.******@***.***> / 979-458-9845 / RDMC 212 On Jun 1, 2023, at 4:35 PM, Demkowicz, Michael J ***@***.******@***.***>> wrote: I’m not familiar with the FIRE algorithm. However, overdamped dynamics is very close to constant increment steepest descent. Glancing over the paper, I see it has standard rigid body dynamics with the only apparent novelty being in the integration scheme. This makes me think it may be possible to do a quick kluge within lammps script where you run 1 tilmestep at a time in a loop. On each step, you set the COM velocity to the COM force and the angular velocity to the torque for each molecule. Then you adjust the tilmestep to limit the max displacement and rotation. This may be a relatively simple thing to try first with no need to alter the source code. — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.******@***.***> / 979-458-9845 / RDMC 212 On May 31, 2023, at 4:07 PM, Steve Plimpton ***@***.******@***.***>> wrote: This Message Is From an External Sender This message came from outside your organization. You should look at this paper for HOOMD. I don't have a copy at the moment, but it has a section on energy min with rigid body constraints. My recollection is it uses FIRE (damped dynamics) not a CG based method. @ndtrung81<https://urldefense.com/v3/__https://github.com/ndtrung81__;!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfcY-bfDCg$> Trung can comment as he is the lead author. I think this is likely a better match to how LAMMPS would use MPI (or GPU) parallelism. Nguyen, T. D., Phillips, C. L., Anderson, J. A., & Glotzer, S. C. (2011). Rigid body constraints realized in massively-parallel molecular dynamics on graphics processing units. Computer Physics Communications, 182(11), 2307-2313. https://doi.org/10.1016/j.cpc.2011.06.005<https://urldefense.com/v3/__https://doi.org/10.1016/j.cpc.2011.06.005__;!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfcpxsTSGI$> — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/lammps/lammps/issues/3801*issuecomment-1570956487__;Iw!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfcZqjwuos$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AD3OZLUS7GDZVE3XVAKK46DXI6XH7ANCNFSM6AAAAAAYVJ665Y__;!!KwNVnqRv!CUqmCFWEN700XuIkdt53HJwi1-N1WNo4DKtPVkPmW0y9aXA4aXisWOJwrUYNPrk_sxh6Ioow5RT9VGW0iDfc3Dy3nWE$>. You are receiving this because you authored the thread.

[athomps]

@demkowicz A lot of the difficulties with minimization in LAMMPS is that atoms may migrate between processors after every energy evaluation. This requires carefully updating all data structures, including those for rigid constraints. It is perfectly doable, but requires some effort.

[demkowicz]

Gotcha. But doesn’t this come up with regular dynamics, too? — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.*** / 979-458-9845 / RDMC 212 On Jun 1, 2023, at 4:50 PM, Aidan Thompson ***@***.***> wrote:  @demkowicz A lot of the difficulties with minimization in LAMMPS is that atoms may migrate between processors after every energy evaluation. This requires carefully updating all data structures, including those for rigid constraints. It is ZjQcmQRYFpfptBannerStart This Message Is From an External Sender This message came from outside your organization. ZjQcmQRYFpfptBannerEnd @demkowicz<https://urldefense.com/v3/__https://github.com/demkowicz__;!!KwNVnqRv!HwHowefE5yRQTv2rVsi4tqSrFQTW3T_SYUyJm0Q7_d88GMKq6G6Uz7KH2taEsV0wP6RCgb7t362KabPoY6_5oTtHSjE$> A lot of the difficulties with minimization in LAMMPS is that atoms may migrate between processors after every energy evaluation. This requires carefully updating all data structures, including those for rigid constraints. It is perfectly doable, but requires some effort. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/lammps/lammps/issues/3801*issuecomment-1572835386__;Iw!!KwNVnqRv!HwHowefE5yRQTv2rVsi4tqSrFQTW3T_SYUyJm0Q7_d88GMKq6G6Uz7KH2taEsV0wP6RCgb7t362KabPoY6_5s3wQhrw$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AD3OZLQJEFWNEJKOMHYMVXLXJEFCBANCNFSM6AAAAAAYVJ665Y__;!!KwNVnqRv!HwHowefE5yRQTv2rVsi4tqSrFQTW3T_SYUyJm0Q7_d88GMKq6G6Uz7KH2taEsV0wP6RCgb7t362KabPoY6_53cnyYlY$>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[athomps]

Yes, but not every timestep, because the halo region is buffered. Also, there is something else. The line search methods, like CG, do backtracking, so they start out with a big step and then shrink it.

[demkowicz]

Yeah, I figured it would be something like that. Thanks for clarifying! — Michael J Demkowicz, Professor Texas A&M Materials Science and Engineering ***@***.******@***.***> / 979-458-9845 / RDMC 212 On Jun 1, 2023, at 4:57 PM, Aidan Thompson ***@***.******@***.***>> wrote: This Message Is From an External Sender This message came from outside your organization. Yes, but not every timestep, because the halo region is buffered. Also, there is something else. The line search methods, like CG, do backtracking, so they start out with a big step and then shrink it. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/lammps/lammps/issues/3801*issuecomment-1572841207__;Iw!!KwNVnqRv!CEj_TGk51mx7edQCGW36sYGda0aKZJiLe2p17mbo8sUMShdcsbKal71bfzl9gu5_HvXcIg02BpTzoXT2FFpcZLx2jjg$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AD3OZLVGE3KCC7OKLYPXUTTXJEF3HANCNFSM6AAAAAAYVJ665Y__;!!KwNVnqRv!CEj_TGk51mx7edQCGW36sYGda0aKZJiLe2p17mbo8sUMShdcsbKal71bfzl9gu5_HvXcIg02BpTzoXT2FFpcrcycYAQ$>. You are receiving this because you were mentioned.

[demkowicz]

Well, I wrote my little poor man's rigid body relaxation algorithm. It works! It takes about 70 iterations (about 1s on a single CPU) to relax 64 tip4p H2O molecules in an ice Ih structure to within about 0.01eV. Script available upon request.

[mm5986]

Well, I wrote my little poor man's rigid body relaxation algorithm. It works! It takes about 70 iterations (about 1s on a single CPU) to relax 64 tip4p H2O molecules in an ice Ih structure to within about 0.01eV. Script available upon request.

This is great! I've been looking for a rigid body relaxation algorithm. Could you please send/post your script.

[demkowicz]

See attached. Note that I have not done anything more than a proof of concept. In particular, I did not do any validation or optimization, so I give absolutely no guarantees on performance. There is also clear room for improvement, e.g., replacing the energy-based convergence criterion with a force/torque-based one.

in.lammps_relax.txt

[mm5986]

Thanks a ton. This will be so helpful.

…

On Thu, Jul 6, 2023, 2:22 AM demkowicz ***@***.***> wrote: See attached. Note that I have not done anything more than a proof of concept. In particular, I did not do any validation or optimization, so I give absolutely no guarantees on performance. There is also clear room for improvement, e.g., replacing the energy-based convergence criterion with a force/torque-based one. in.lammps_relax.txt <https://github.com/lammps/lammps/files/11964976/in.lammps_relax.txt> — Reply to this email directly, view it on GitHub <#3801 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BBCCWGABLH2AYGTZOVOBP3TXOZKQ5ANCNFSM6AAAAAAYVJ665Y> . You are receiving this because you commented.Message ID: ***@***.***>