Issue with writing virtual sites via a lone pair

[CalCraven]

For forcefields such as tip4p, we can read the forcefield, and virtual sites are treated just the same as a site. The virtual site is then assigned an atomtype just the same via the forcefield, and this is great. However, for our writers, these virtual sites need to be identified (and differentiated from coarse grained systems) in order to be handled specially. For example, Gromacs has a [ dummies ] section that needs to be written if a lone pair is included. Currently, since we have no special implementation for it, the top writer will throw an error.

We need to add some logic into the writers that can handle virtual sites, potentially some changes logic into the gmso.Site to assist with that, and make sure any readers that could have virtual sites would also properly bring that into a gmso.Topology.

import mbuild
from gmso.parameterization import apply

OPC_XML = 'OPC4_GMSO.xml'

opc_wat = mb.load('OPC.mol2')

N_COMPOUNDS=100
X_LENGTH = 5.15
startingBox=mb.fill_box(compound=opc_wat ,n_compounds=2,box=[5,5,5])
top = startingBox.to_gmso()

gmso_ff = gmso.ForceField('OPC4_GMSO.xml')
parameterized_top = apply(top, gmso_ff, identify_connections=True, remove_untyped=True) 

parameterized_top.save('init_gmso.gro', overwrite=True)
parameterized_top.save('init_gmso.top', overwrite=True)

Errors as:

UnexpectedCharacters                      Traceback (most recent call last)
File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/lexer.py:486, in ContextualLexer.lex(self, lexer_state, parser_state)
    485 last_token = lexer_state.last_token  # Save last_token. Calling root_lexer.next_token will change this to the wrong token
--> 486 token = self.root_lexer.next_token(lexer_state, parser_state)
    487 raise UnexpectedToken(token, e.allowed, state=parser_state, token_history=[last_token], terminals_by_name=self.root_lexer.terminals_by_name)

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/lexer.py:398, in TraditionalLexer.next_token(self, lex_state, parser_state)
    397         allowed = {"<END-OF-FILE>"}
--> 398     raise UnexpectedCharacters(lex_state.text, line_ctr.char_pos, line_ctr.line, line_ctr.column,
    399                                allowed=allowed, token_history=lex_state.last_token and [lex_state.last_token],
    400                                state=parser_state, terminals_by_name=self.terminals_by_name)
    402 value, type_ = res

UnexpectedCharacters: No terminal matches '_' in the current parser context, at line 1 col 2

[_LP;X1][O;X3](H)(H)
 ^
Expected one of: 
	* BANG
	* SEMICOLON
	* LABEL
	* "R"
	* "$("
	* COMMA
	* LPAR
	* X
	* RPAR
	* R
	* LSQB
	* RSQB
	* NUM
	* STAR
	* AMPERSAND
	* SYMBOL
	* HASH

Previous tokens: Token('LSQB', '[')


During handling of the above exception, another exception occurred:

UnexpectedCharacters                      Traceback (most recent call last)
Cell In[23], line 46
     44 print(f"Found {len(parameterized_top.atom_types)} in topology")
     45 parameterized_top.save('init_gmso.gro', overwrite=True)
---> 46 parameterized_top.save('init_gmso.top', overwrite=True)

File ~/Dropbox/Mac/Documents/Vanderbilt/Research/MoSDeF/gmso/gmso/core/topology.py:1490, in Topology.save(self, filename, overwrite, **kwargs)
   1487 from gmso.formats import SaversRegistry
   1489 saver = SaversRegistry.get_callable(filename.suffix)
-> 1490 saver(self, filename, **kwargs)

File ~/Dropbox/Mac/Documents/Vanderbilt/Research/MoSDeF/gmso/gmso/formats/top.py:102, in write_top(top, filename, top_vars)
     81 out_file.write(
     82     "[ atomtypes ]\n"
     83     "; name\t"
   (...)
     89     "epsilon\n"
     90 )
     92 for atom_type in top.atom_types(PotentialFilters.UNIQUE_NAME_CLASS):
     93     out_file.write(
     94         "{0:12s}"
     95         "{1:4s}"
     96         "{2:12.5f}"
     97         "{3:12.5f}\t"
     98         "{4:4s}"
     99         "{5:12.5f}"
    100         "{6:12.5f}\n".format(
    101             atom_type.name,
--> 102             str(_lookup_atomic_number(atom_type)),
    103             atom_type.mass.in_units(u.amu).value,
    104             atom_type.charge.in_units(u.elementary_charge).value,
    105             "A",
    106             atom_type.parameters["sigma"].in_units(u.nanometer).value,
    107             atom_type.parameters["epsilon"]
    108             .in_units(u.Unit("kJ/mol"))
    109             .value,
    110         )
    111     )
    113 # Define unique molecule by name only
    114 unique_molecules = _get_unique_molecules(top)

File ~/Dropbox/Mac/Documents/Vanderbilt/Research/MoSDeF/gmso/gmso/formats/top.py:381, in _lookup_atomic_number(atom_type)
    379 """Look up an atomic_number based on atom type information, 0 if non-element type."""
    380 try:
--> 381     element = element_by_atom_type(atom_type)
    382     return element.atomic_number
    383 except GMSOError:

File ~/Dropbox/Mac/Documents/Vanderbilt/Research/MoSDeF/gmso/gmso/core/element.py:313, in element_by_atom_type(atom_type, verbose)
    311     matched_element = element_by_symbol(atom_type.name, verbose=verbose)
    312 if matched_element is None and atom_type.definition:
--> 313     matched_element = element_by_smarts_string(
    314         atom_type.definition, verbose=verbose
    315     )
    317 if matched_element is None:
    318     raise GMSOError(
    319         f"Failed to find an element from atom type"
    320         "{atom_type} with "
   (...)
    323         "definition: {atom_type.definition}"
    324     )

File ~/Dropbox/Mac/Documents/Vanderbilt/Research/MoSDeF/gmso/gmso/core/element.py:260, in element_by_smarts_string(smarts_string, verbose)
    256 SMARTS = foyer.smarts.SMARTS
    258 PARSER = SMARTS()
--> 260 symbols = PARSER.parse(smarts_string).iter_subtrees_topdown()
    262 first_symbol = None
    263 for symbol in symbols:

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/foyer/smarts.py:86, in SMARTS.parse(self, smarts_string)
     84 def parse(self, smarts_string):
     85     """Convert SMARTS string to parsed grammar object."""
---> 86     return self.PARSER.parse(smarts_string)

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/lark.py:581, in Lark.parse(self, text, start, on_error)
    563 def parse(self, text, start=None, on_error=None):
    564     """Parse the given text, according to the options provided.
    565 
    566     Parameters:
   (...)
    579 
    580     """
--> 581     return self.parser.parse(text, start=start, on_error=on_error)

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/parser_frontends.py:106, in ParsingFrontend.parse(self, text, start, on_error)
    104 stream = text if self.skip_lexer else LexerThread(self.lexer, text)
    105 kw = {} if on_error is None else {'on_error': on_error}
--> 106 return self.parser.parse(stream, chosen_start, **kw)

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/parsers/lalr_parser.py:41, in LALR_Parser.parse(self, lexer, start, on_error)
     39 def parse(self, lexer, start, on_error=None):
     40     try:
---> 41         return self.parser.parse(lexer, start)
     42     except UnexpectedInput as e:
     43         if on_error is None:

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/parsers/lalr_parser.py:171, in _Parser.parse(self, lexer, start, value_stack, state_stack, start_interactive)
    169 if start_interactive:
    170     return InteractiveParser(self, parser_state, parser_state.lexer)
--> 171 return self.parse_from_state(parser_state)

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/parsers/lalr_parser.py:188, in _Parser.parse_from_state(self, state)
    186     except NameError:
    187         pass
--> 188     raise e
    189 except Exception as e:
    190     if self.debug:

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/parsers/lalr_parser.py:178, in _Parser.parse_from_state(self, state)
    176 try:
    177     token = None
--> 178     for token in state.lexer.lex(state):
    179         state.feed_token(token)
    181     end_token = Token.new_borrow_pos('$END', '', token) if token else Token('$END', '', 0, 1, 1)

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/lexer.py:489, in ContextualLexer.lex(self, lexer_state, parser_state)
    487     raise UnexpectedToken(token, e.allowed, state=parser_state, token_history=[last_token], terminals_by_name=self.root_lexer.terminals_by_name)
    488 except UnexpectedCharacters:
--> 489     raise e

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/lexer.py:478, in ContextualLexer.lex(self, lexer_state, parser_state)
    476     while True:
    477         lexer = self.lexers[parser_state.position]
--> 478         yield lexer.next_token(lexer_state, parser_state)
    479 except EOFError:
    480     pass

File ~/miniconda3/envs/gmso-dev/lib/python3.11/site-packages/lark/lexer.py:398, in TraditionalLexer.next_token(self, lex_state, parser_state)
    396     if not allowed:
    397         allowed = {"<END-OF-FILE>"}
--> 398     raise UnexpectedCharacters(lex_state.text, line_ctr.char_pos, line_ctr.line, line_ctr.column,
    399                                allowed=allowed, token_history=lex_state.last_token and [lex_state.last_token],
    400                                state=parser_state, terminals_by_name=self.terminals_by_name)
    402 value, type_ = res
    404 if type_ not in self.ignore_types:

UnexpectedCharacters: No terminal matches '_' in the current parser context, at line 1 col 2

[_LP;X1][O;X3](H)(H)
 ^
Expected one of: 
	* BANG
	* R
	* LABEL
	* "R"
	* STAR
	* "$("
	* X
	* SYMBOL
	* HASH

Previous tokens: Token('LSQB', '[')

[CalCraven]

#196 Is also discussing an implementation long ago

[CalCraven]

We also have a few PR's that have been opened previously to look at this:
VirtualSite class: #197
Beads, VirtualSites, DummyAtoms: #523

[CalCraven]

GROMACS implementation is here, OPENFF implementation is here, LAMMPS implementation here