search

mosdef-hub / mbuild

A hierarchical, component based molecule builder
https://mbuild.mosdef.org
Other
171 stars 80 forks source link

FoyerError: Found no types for atom 0 (None) when saving LAMMPS file #1101

Closed JiaYuanChng closed 1 year ago

JiaYuanChng commented 1 year ago

Hello,

I have not been able to fix an issue associated with saving a lammps file using the .save function.

My code for constructing united-atom type alkane molecules, and the forcefield .xml file are attached below. While everything worked fine for ethane (ie. chain_length=1) and lammps file was successfully created, it stopped working when chain_length=2. Visualization of propane looked fine, inspection using test_propane.labels also seems to look fine.

Have I overlooked something here?

Thank you, JY

import mbuild as mb
import warnings
warnings.filterwarnings('ignore')

class CH3(mb.Compound):
    def __init__(self):
        super(CH3, self).__init__()

        self.add(mb.Particle(name='_CH3', pos=[0,0,0]), label='CH3')
        self.add(mb.Port(anchor=self['CH3'][0]), label='up')
        self['up'].translate([0.154/2, 0, 0])

class CH2(mb.Compound):
    def __init__(self):
        super(CH2, self).__init__()

        self.add(mb.Particle(name='_CH2', pos=[0,0,0]), label='CH2')
        self.add(mb.Port(anchor=self['CH2'][0]), label='up')
        self.add(mb.Port(anchor=self['CH2'][0]), label='down')
        self['up'].translate([0.154/2, 0, 0])
        self['down'].translate([-0.154/2, 0, 0])

class Alkane(mb.Compound):
    def __init__(self, chain_length=1):
        super(Alkane, self).__init__()

        last_monomer = CH3()
        self.add(last_monomer)

        # loop over n = chain_length CH2's to add to Alkane
        for _ in range(chain_length-1):
            current_monomer = CH2()
            mb.force_overlap(move_this=current_monomer,
                             from_positions=current_monomer['down'],
                             to_positions=last_monomer['up'])
            self.add(current_monomer)
            last_monomer=current_monomer

        last_cap = CH3()
        mb.force_overlap(move_this=last_cap,
                         from_positions=last_cap['up'],
                         to_positions=last_monomer['up'])
        self.add(last_cap)

Ethane (works fine):

test_ethane = Alkane(chain_length=1)
test_ethane.save('UA_test_ethane.lammps', forcefield_files='./trappe.xml', overwrite=True)

Propane:

test_propane = Alkane(chain_length=2)
test_propane.save('UA_test_propane.lammps', forcefield_files='./trappe.xml', overwrite=True)

Error:

---------------------------------------------------------------------------
FoyerError                                Traceback (most recent call last)
Cell In[53], line 1
----> 1 test_propane.save('UA_test_propane.lammps', forcefield_files='./UA_ff/trappe.xml', overwrite=True)

File ~/.conda/envs/mosdef/lib/python3.9/site-packages/mbuild/compound.py:2683, in Compound.save(self, filename, show_ports, forcefield_name, forcefield_files, forcefield_debug, box, overwrite, residues, combining_rule, foyer_kwargs, parmed_kwargs, **kwargs)
   2576 def save(
   2577     self,
   2578     filename,
   (...)
   2589     **kwargs,
   2590 ):
   2591     """Save the Compound to a file.
   2592 
   2593     Parameters
   (...)
   2681     formats.json_formats.compound_to_json : Write to a json file
   2682     """
-> 2683     conversion.save(
   2684         self,
   2685         filename,
   2686         show_ports,
   2687         forcefield_name,
   2688         forcefield_files,
   2689         forcefield_debug,
   2690         box,
   2691         overwrite,
   2692         residues,
   2693         combining_rule,
   2694         foyer_kwargs,
   2695         parmed_kwargs,
   2696         **kwargs,
   2697     )

File ~/.conda/envs/mosdef/lib/python3.9/site-packages/mbuild/conversion.py:1072, in save(compound, filename, show_ports, forcefield_name, forcefield_files, forcefield_debug, box, overwrite, residues, combining_rule, foyer_kwargs, parmed_kwargs, **kwargs)
   1070 if not foyer_kwargs:
   1071     foyer_kwargs = {}
-> 1072 structure = ff.apply(structure, **foyer_kwargs)
   1073 if structure.combining_rule != combining_rule:
   1074     warn(
   1075         f"Overwriting forcefield-specified combining rule ({combining_rule})"
   1076         f"to new combining rule ({combining_rule})."
   (...)
   1079         "Consider directly changing the metadata of the Forcefield."
   1080     )

File ~/.conda/envs/mosdef/lib/python3.9/site-packages/foyer/forcefield.py:780, in Forcefield.apply(self, structure, references_file, use_residue_map, assert_bond_params, assert_angle_params, assert_dihedral_params, assert_improper_params, verbose, *args, **kwargs)
    777     if isinstance(structure, mb.Compound):
    778         structure = structure.to_parmed(**kwargs)
--> 780 typemap = self.run_atomtyping(
    781     structure, use_residue_map=use_residue_map, **kwargs
    782 )
    784 self._apply_typemap(structure, typemap)
    786 structure = self.parametrize_system(
    787     structure=structure,
    788     references_file=references_file,
   (...)
    795     **kwargs,
    796 )

File ~/.conda/envs/mosdef/lib/python3.9/site-packages/foyer/forcefield.py:853, in Forcefield.run_atomtyping(self, structure, use_residue_map, **kwargs)
    848     if (
    849         structure.residues[res_id].name
    850         not in residue_map.keys()
    851     ):
    852         tmp_res = _structure_from_residue(res, structure)
--> 853         typemap = find_atomtypes(tmp_res, forcefield=self)
    854         residue_map[res.name] = typemap
    856 typemap = _unwrap_typemap(structure, residue_map)

File ~/.conda/envs/mosdef/lib/python3.9/site-packages/foyer/atomtyper.py:145, in find_atomtypes(structure, forcefield, max_iter)
    142 rules = subrules
    144 _iterate_rules(rules, topology_graph, typemap, max_iter=max_iter)
--> 145 _resolve_atomtypes(topology_graph, typemap)
    147 return typemap

File ~/.conda/envs/mosdef/lib/python3.9/site-packages/foyer/atomtyper.py:225, in _resolve_atomtypes(topology_graph, typemap)
    218     raise FoyerError(
    219         "Found multiple types for atom {} ({}): {}.".format(
    220             atom_id, atoms[atom_id].atomic_number, atomtype
    221         )
    222     )
    223 else:
--> 225     raise FoyerError(
    226         "Found no types for atom {} ({}).".format(
    227             atom_id, atoms[atom_id].atomic_number
    228         )
    229     )

FoyerError: Found no types for atom 0 (None).

Force field .xml file:

<ForceField>
  <AtomTypes>
    <Type name="CH4" class="Csp3" element="_CH4" mass="16.04"
          def="[_CH4;X0]" desc="UA CH4"/>
    <Type name="CH3_alk" class="Csp3" element="_CH3" mass="15.02348"
          def="[_CH3;X1](_CH3)" desc="UA CH3 alkane"/>
    <Type name="CH2_alk" class="Csp3" element="_CH2" mass="14.027"
          def="[_CH2;X2]([_CH3,_CH2])[_CH3,_CH2]" desc="UA CH2 alkane"/>
  </AtomTypes>
  <HarmonicBondForce>
      <Bond type1="CH3_alk" type2="CH3_alk" length="0.154" k="5000"/>
      <Bond type1="CH3_alk" type2="CH2_alk" length="0.154" k="5000"/>
  </HarmonicBondForce>
  <HarmonicAngleForce>
      <Angle type1="CH3_alk" type2="CH2_alk" type3="CH3_alk" angle="1.989675" k="519.651"/>
      <Angle type1="CH3_alk" type2="CH2_alk" type3="CH2_alk" angle="1.989675" k="519.651"/>
      <Angle type1="CH2_alk" type2="CH2_alk" type3="CH2_alk" angle="1.989675" k="519.651"/>
  </HarmonicAngleForce>
  <RBTorsionForce>
      <Proper type1="CH3_alk" type2="CH2_alk" type3="CH2_alk" type4="CH2_alk"
        c0="2.9288" c1="-1.4644" c2="0.2092" c3="-1.6736" c4="0.0" c5="0.0"/>
      <Proper type1="CH3_alk" type2="CH2_alk" type3="CH2_alk" type4="CH3_alk"
        c0="2.9288" c1="-1.4644" c2="0.2092" c3="-1.6736" c4="0.0" c5="0.0"/>
      <Proper type1="CH2_alk" type2="CH2_alk" type3="CH2_alk" type4="CH2_alk"
        c0="2.9288" c1="-1.4644" c2="0.2092" c3="-1.6736" c4="0.0" c5="0.0"/>
  </RBTorsionForce>
  <NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
    <Atom type="CH4" charge="0.00" sigma="0.373" epsilon="1.231"/>
    <Atom type="CH3_alk" charge="0.00" sigma="0.375" epsilon="0.815"/>
    <Atom type="CH2_alk" charge="0.00" sigma="0.395" epsilon="0.382"/>
  </NonbondedForce>
</ForceField>
JiaYuanChng commented 1 year ago

Blunder resolved: Inaccurate description for CH3_alk: def="[_CH3;X1]([_CH3,_CH2])"