Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support for non-element for GROMACS writer #817

Open
daico007 opened this issue Apr 12, 2024 · 5 comments
Open

Support for non-element for GROMACS writer #817

daico007 opened this issue Apr 12, 2024 · 5 comments
Labels
bug Something isn't working core feature enhancement New feature or request writer writer related

Comments

@daico007
Copy link
Member

Current GROMACS writer didn't take into the case of system with non-element site. This leads to a wide range of errors that may not be obvious to user. The method need to be revised to take those into account + add more unit tests.

@bc118
Copy link
Contributor

bc118 commented Apr 12, 2024

I found the line that is causing the issues.
In the "gmso location/gmso/parameterization

The atomic number needs to be zero (0) for all dummy atoms. For example, for TIP4P, that is what parmed uses to set the water settle/rigid, and if the dummy atom is not there, it will not set settle/fix rigid correctly. Simply loading the gmso printed .top file in parmed and rewriting it back out with parmed fixes the TIP4P to have the correct setup.

Parmed rewriting for TIP4P, is a good place to see how the it should be setup as it does what we want, but we also need to exclude things that the parmed writer does. This is how you can see what it does:

from parmed.gromacs import GromacsTopologyFile

# load into parmed and rewrite to fix tip4op top file
parmed_top = GromacsTopologyFile('init.top')
print(parmed_top.topology)
#parmed_top.createSystem( rigidWater=True)
parmed_top.save('init.top', overwrite=True)

The final gromacs output for TIP4P will be attached for both the 'flexible and 'rigid' TIP4P model (see attachments in this issue):

I am not sure what all tests you want to include so I will leave it at this, unless you want me to push a PR

This fixes the issue with not loading Dummy atoms or _CH3 (i.e., beads) to be written in the GMSO GROMACS .top writer:
In the "gmso location/gmso/parameterization

def _lookup_atomic_number(atom_type):
    """Look up an atomic_number based on atom type information, 0 if non-element type."""
    try:
        element = element_by_atom_type(atom_type)
        return element.atomic_number
    except:
        return 0


def _lookup_element_symbol(atom_type):
    """Look up an atomic_number based on atom type information, 0 if non-element type."""
    try:
        element = element_by_atom_type(atom_type)
        return element.symbol
    except:
        return "X"

was this (basically remove the "GMSOError" from them:

def _lookup_atomic_number(atom_type):
    """Look up an atomic_number based on atom type information, 0 if non-element type."""
    try:
        element = element_by_atom_type(atom_type)
        return element.atomic_number
    except GMSOError:
        return 0


def _lookup_element_symbol(atom_type):
    """Look up an atomic_number based on atom type information, 0 if non-element type."""
    try:
        element = element_by_atom_type(atom_type)
        return element.symbol
    except GMSOError:
        return "X"

@bc118 bc118 added bug Something isn't working enhancement New feature or request core feature writer writer related labels Apr 18, 2024
@dyukovsm
Copy link

dyukovsm commented Apr 18, 2024

I am attaching a validated tip4p-ew topology here for both the 'flexible' and 'rigid' TIP4P models. From talking with @dyukovsm , the GROMACS manual says to include or do it a different way. However, @dyukovsm has tested this with Jeff and found that the below produces the correct results, not what the GROMACS manual does. @dyukovsm is creating and issue with GROMACS on this topic

The 'rigid' TIP4P :

;
;   File init.top  was generated
;   By user: brad (1000)
;   On host:z840-1-brad 
;   At date:Sat. April  6 15:45:07 2024 
;
;   This is a standalone topology file
;
;   Created by:
;   ParmEd:       VERSION4.2.2 
;   Executable:   project.py
;   Library dir: xxxx
;   Command line:
;     xxxxx
;

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               3               yes             0.5          0.5     

[ atomtypes ]
; name    at.num    mass    charge ptype  sigma      epsilon
OW             8  15.999400  0.00000000  A        0.31644        0.68095
HW             1   1.008000  0.52422000  A              0              0
_LP            0   0.000000  -1.04844000  A              0              0


[ moleculetype ]
; Name            nrexcl
TIP4          3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr    charge       mass  typeB    chargeB      massB
; residue    1 TIP4 rtp TIP4 q 0.0
    1         OW      1   TIP4      O      1 0.00000000  15.999400   ; qtot 0.000000
    2         HW      1   TIP4      H      2 0.52422000   1.008000   ; qtot 0.524220
    3         HW      1   TIP4      H      3 0.52422000   1.008000   ; qtot 1.048440
    4        _LP      1   TIP4    _LP      4 -1.04844000   0.000000   ; qtot 0.000000

[ bonds ]
;    ai     aj funct         c0         c1         c2         c3
      1      3     1   0.09572 502416.000000
      1      2     1   0.09572 502416.000000

[ angles ]
;    ai     aj     ak funct         c0         c1         c2         c3
      2      1      3     1   104.5200000 628.020000
      2      1      4     1   52.2600000 418.680000
      3      1      4     1   52.2600000 418.680000

[ settles ]
; i     funct   doh     dhh
1     1   0.09572000   0.15139007


[ virtual_sites3 ]
; Site  from                   funct
4     1    3    2    1      0.106677  0.106677

[ exclusions ]
1  2  3  4
2  1  4  3  4
3  1  4  2  4
4  1  2  3

[ system ]
; Name
Generic title

[ molecules ]
; Compound       #mols
WAT               4500

'flexible' TIP4P = Removing the below section from the 'rigid' TIP4P:

[ settles ]
; i     funct   doh     dhh
1     1   0.09572000   0.15139007

@bc118
Copy link
Contributor

bc118 commented Apr 18, 2024

Attached is a mol2 and GMSO FF to recreate and test this

tip4p_ew_mol2_and_ff.zip

@dyukovsm
Copy link

python script write_gmx_tip4p.py reads TIP4.mol2 and using foyer reads TIP4P-2005.xml, and writes the UNCONSTRAINED_init.top file. It's pretty dirty and didn't have options for constraints (that I know of), so I read the UNCONSTRAINED_init.top with parmed, save it suing parmed as IFNDEF_init.top, and copy the file using native python text editor without ifndef stuff, because those seem to not run correctly.

TEST.zip

@CalCraven
Copy link
Contributor

Thanks for all this info @bc118 @dyukovsm! This is really helpful having specific files and examples to work with when writing unit tests. I would say because the various water methods are so ubiquitous and requested for, let's make some changes in #771 for handling this case specifically, instead of a general method for non-element rigid body virtual sites.

On that note, there are two things we need to be able to handle in that PR.

  1. rigid water. Since we don't have a way to specify this in the xml, let's just make it an argument to the gromacs writer.
    rigid_water=True.
    This will write the settles section, and is easily toggleable. Can leave in the current logic checking for HOH molecules, but remove the rigid check for a molecule label until we can rework that entire class, as discussed in that PR.
  2. Virtual site water. This is even more tricky, but also extremely important. Again, a simple gromacs specific solution is to pass in an argument to the writer:
    3_site_virtual_types = ["_LP"]
    That way a user can specify what needs to be looked for, and easily identify when to write out the virtual sites section. We will only be able to handle the 3 site type, but that may be as good a bandaid for now as we can expect until a more general solution that is flexible to different writers.

Again, these solutions are bandaids. Ideally this information is entirely contained in the xml file. So we'll need methods to extract a virtual site, calculate the general parameters, and port that information to different writers. Again doable, but will need an overhaul of the xml formatting that can happen at a later date.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working core feature enhancement New feature or request writer writer related
Projects
None yet
Development

No branches or pull requests

4 participants