Incomplete accounting for fragment energy in CBMC moves

[rwsmith7531]

Wouldn't it be possible for a non-ring fragment to have improper torsions due to non-ring sp2 atoms, whether represented as harmonic impropers or improper dihedrals? If so, the code very much looks like they wouldn't be taken into account for the dE_frag in CBMC moves, since the fragment growth routines only check the fragment energy for ring fragments, apparently making the assumption that only bond angle energy influenced the fragment conformation for non-ring fragments. This assumption is probably typically correct, but not necessarily always correct. The Atom_Displacement subroutine called by the fragment driver computes not only angle energy, but also dihedral energy (which may include an improper torsion input as a dihedral), improper energy, and nonbonded energy. Even fragments without any sp2 hybridization could have nonzero nonbonded intramolecular energy if the 1-2 or 1-3 intramolecular scaling factors are nonzero, which may be rare but is allowed by Cassandra.

I also noticed that all angle energy is included in both dE_frag and dE (via dE_intra) for CBMC moves (and Widom insertions), and the angle energy would just cancel out when computing the move probability or widom_var as dE_frag is subtracted from dE. I don't think angle energy should be computed at all for those moves because it doesn't matter.

[rwsmith7531]

I don't think any of the test suite cases include impropers either, except for the test specifically for impropers and the trajectory reader test with BMIM-PF6 (I think there's an improper input as a dihedral in BMIM), so none of the tests that include them perform CBMC moves on species that have them, so this wouldn't be caught by the test suite.

[emarinri]

Hello Ryan,

I think you raise very good points. The key assumption of the fragment-based sampling is that the fragment configuration in the gas phase (i.e., the angles) is the same as the fragment configuration in the liquid phase, since the angles are stiff degrees of freedom and thus are essentially insensitive to the environment. We currently assume that the angle energy arises only due to the angle harmonic potential.

As it is typical in organic fluids, non-bonded 1-2 or 1-3 interactions are 0.0. I assume in more "exotic" force fields (at least exotic relative to our typical use cases) these could be non-zero. For these cases, could angles still be considered stiff? or would the assumption of separation of stiff and soft degrees of freedom (dihedrals) would break? My first guess is that our basic assumption should hold. VdW energies for 1-2 and 1-3 interactions should be very strong and lead to even stiffer angles due to short distances. If the basic assumption still holds, I think the total fragment energy should include non-bonded 1-2 and 1-3 as well, and such energy should also go in the acceptance rule.

Perhaps someone is aware of cases where 1-2 and 1-3 interactions are not zero for organic force fields? what are examples of 1-2 and 1-3 non-bonded interactions that are not zero in general? (some special force fields for biomolecules? materials? coarse-grained models? custom force fields?)

For the case of improper torsional angles whose constituent atoms are the same as those that make up "typical" fragment (e.g. a central atom connected to two or three atoms), this improper energy should be included in the fragment energy as well. Again, the acceptance rule should include this energy. My guess is that the separation of stiff and soft degrees of freedom would still hold.

Cassandra/Src/move_atom.f90

Lines 253 to 258 in aeefdf1

	CALL Compute_Molecule_Angle_Energy(lm,is,E_angle_o)
	CALL Compute_Molecule_Dihedral_Energy(lm,is,E_dihed_o)
	CALL Compute_Molecule_Improper_Energy(lm,is,E_improper_o)
	dE_intra = dE_intra - E_angle_o - E_dihed_o - E_improper_o
	CALL Compute_Atom_Nonbond_Energy(iatom,lm,is,&

As you mention, the Atom_Displacement move computes dihedral, improper and non-bonded energies. This move is just another move in the available moves the code provides. It could be used in a typical simulation, in which you might need to compute such energies that would go in the acceptance rule. However, we typically use it just in the fragment library generation, as it is part of the usual workflow.

I don't think any of the test suite cases include impropers either, except for the test specifically for impropers and the trajectory reader test with BMIM-PF6 (I think there's an improper input as a dihedral in BMIM), so none of the tests that include them perform CBMC moves on species that have them, so this wouldn't be caught by the test suite.

We do not. I am not sure if we test fragment energies in any way.

[rwsmith7531]

I also realized that when ring fragment energy is computed with del_flag = .TRUE. (i.e. deletion moves, part of regrowth, part of mol swap), the electrostatic portion of the energy it computes won't be the same as the electrostatic energy used to generate the ring conformation during the ring fragment library generation simulation if the charge sum style is anything other than simple minimum image, because the fragment library generation simulations use the minimum image charge sum style (that's how the fragment library generation python script writes the fragment library generation simulation input files). This problem doesn't occur when del_flag is .FALSE. because the intrafragment ring fragment energy is taken from the fragment library. I think this problem could be fixed by setting igas_flag = .TRUE. at the beginning of Compute_Ring_Fragment_Energy. igas_flag currently has no way of being set to .TRUE. anywhere in Cassandra, and only appears where it's declared, where it's initialized to .FALSE., and when it's checked by the electrostatic portion of Compute_AtomPair_Energy (it would cause it to behave as if the charge sum style were simple cutoff if it were .TRUE., but it's never .TRUE.).

…

On Fri, Sep 29, 2023, 6:55 PM Eliseo Marin-Rimoldi ***@***.***> wrote: Hello Ryan, I think you raise very good points. The key assumption of the fragment-based sampling is that the fragment configuration in the gas phase (i.e., the angles) is the same as the fragment configuration in the liquid phase, since the angles are stiff degrees of freedom and thus are essentially insensitive to the environment. We currently assume that the angle energy arises only due to the angle harmonic potential. As it is typical in organic fluids, non-bonded 1-2 or 1-3 interactions are 0.0. I assume in more "exotic" force fields (at least exotic relative to our typical use cases) these could be non-zero. For these cases, could angles still be considered stiff? or would the assumption of separation of stiff and soft degrees of freedom (dihedrals) would break? My first guess is that our basic assumption should hold. VdW energies for 1-2 and 1-3 interactions should be very strong and lead to even stiffer angles due to short distances. If the basic assumption still holds, I think the total fragment energy should include non-bonded 1-2 and 1-3 as well, and such energy should also go in the acceptance rule. Perhaps someone is aware of cases where 1-2 and 1-3 interactions are not zero for organic force fields? what are examples of 1-2 and 1-3 non-bonded interactions that are not zero? (some special force fields for biomolecules? materials? custom force fields?) For the case of improper torsional angles whose constituent atoms are the same as those that make up "typical" fragment (e.g. a central atom connected to two or three atoms), this improper energy should be included in the fragment energy as well. Again, the acceptance rule should include this energy. My guess is that the separation of stiff and soft degrees of freedom would still hold. So, I guess these should probably be included in the acceptance rules of the fragment regrowth moves, in case someone ever needs to simulate such systems. Most of the times, they should equal to zero. https://github.com/MaginnGroup/Cassandra/blob/aeefdf11cc2b80d3541177a56db48e8dcc52b5e0/Src/move_atom.f90#L253-L258 As you mention, the Atom_Displacement move computes dihedral, improper and non-bonded energies. This move is just another move in the available moves the code provides. It could be used in a typical simulation, in which you might need to compute such energies that would go in the acceptance rule. However, we typically use it just in the fragment library generation, as it is part of the usual workflow. I don't think any of the test suite cases include impropers either, except for the test specifically for impropers and the trajectory reader test with BMIM-PF6 (I think there's an improper input as a dihedral in BMIM), so none of the tests that include them perform CBMC moves on species that have them, so this wouldn't be caught by the test suite. We do not. I am not sure if we test fragment energies in any way. — Reply to this email directly, view it on GitHub <#148 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG6B4SPAVYDGPXNK53S35HDX45GVXANCNFSM6AAAAAA5M6ZYMM> . You are receiving this because you were assigned.Message ID: ***@***.***>

[rwsmith7531]

We would also have to make igas_flag override DSF charge sum style, because it currently does not.

…

On Sat, Sep 30, 2023, 10:58 AM Ryan Smith ***@***.***> wrote: I also realized that when ring fragment energy is computed with del_flag = .TRUE. (i.e. deletion moves, part of regrowth, part of mol swap), the electrostatic portion of the energy it computes won't be the same as the electrostatic energy used to generate the ring conformation during the ring fragment library generation simulation if the charge sum style is anything other than simple minimum image, because the fragment library generation simulations use the minimum image charge sum style (that's how the fragment library generation python script writes the fragment library generation simulation input files). This problem doesn't occur when del_flag is .FALSE. because the intrafragment ring fragment energy is taken from the fragment library. I think this problem could be fixed by setting igas_flag = .TRUE. at the beginning of Compute_Ring_Fragment_Energy. igas_flag currently has no way of being set to .TRUE. anywhere in Cassandra, and only appears where it's declared, where it's initialized to .FALSE., and when it's checked by the electrostatic portion of Compute_AtomPair_Energy (it would cause it to behave as if the charge sum style were simple cutoff if it were .TRUE., but it's never .TRUE.). On Fri, Sep 29, 2023, 6:55 PM Eliseo Marin-Rimoldi < ***@***.***> wrote: > Hello Ryan, > > I think you raise very good points. The key assumption of the > fragment-based sampling is that the fragment configuration in the gas phase > (i.e., the angles) is the same as the fragment configuration in the liquid > phase, since the angles are stiff degrees of freedom and thus are > essentially insensitive to the environment. We currently assume that the > angle energy arises only due to the angle harmonic potential. > > As it is typical in organic fluids, non-bonded 1-2 or 1-3 interactions > are 0.0. I assume in more "exotic" force fields (at least exotic relative > to our typical use cases) these could be non-zero. For these cases, could > angles still be considered stiff? or would the assumption of separation of > stiff and soft degrees of freedom (dihedrals) would break? My first guess > is that our basic assumption should hold. VdW energies for 1-2 and 1-3 > interactions should be very strong and lead to even stiffer angles due to > short distances. If the basic assumption still holds, I think the total > fragment energy should include non-bonded 1-2 and 1-3 as well, and such > energy should also go in the acceptance rule. > > Perhaps someone is aware of cases where 1-2 and 1-3 interactions are not > zero for organic force fields? what are examples of 1-2 and 1-3 non-bonded > interactions that are not zero? (some special force fields for > biomolecules? materials? custom force fields?) > > For the case of improper torsional angles whose constituent atoms are the > same as those that make up "typical" fragment (e.g. a central atom > connected to two or three atoms), this improper energy should be included > in the fragment energy as well. Again, the acceptance rule should include > this energy. My guess is that the separation of stiff and soft degrees of > freedom would still hold. > > So, I guess these should probably be included in the acceptance rules of > the fragment regrowth moves, in case someone ever needs to simulate such > systems. Most of the times, they should equal to zero. > > > https://github.com/MaginnGroup/Cassandra/blob/aeefdf11cc2b80d3541177a56db48e8dcc52b5e0/Src/move_atom.f90#L253-L258 > > As you mention, the Atom_Displacement move computes dihedral, improper > and non-bonded energies. This move is just another move in the available > moves the code provides. It could be used in a typical simulation, in which > you might need to compute such energies that would go in the acceptance > rule. However, we typically use it just in the fragment library generation, > as it is part of the usual workflow. > > I don't think any of the test suite cases include impropers either, > except for the test specifically for impropers and the trajectory reader > test with BMIM-PF6 (I think there's an improper input as a dihedral in > BMIM), so none of the tests that include them perform CBMC moves on species > that have them, so this wouldn't be caught by the test suite. > > We do not. I am not sure if we test fragment energies in any way. > > — > Reply to this email directly, view it on GitHub > <#148 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AG6B4SPAVYDGPXNK53S35HDX45GVXANCNFSM6AAAAAA5M6ZYMM> > . > You are receiving this because you were assigned.Message ID: > ***@***.***> >