Lammps lj unit charges (#1086)

* Produce n=2 when ports are not capped

* Add option to not scale charges in lj units, and update docs for consistency

* Update mbuild/formats/lammpsdata.py

Co-authored-by: Co Quach <43968221+daico007@users.noreply.github.com>

---------

Co-authored-by: Co Quach <43968221+daico007@users.noreply.github.com>

Author: jaclark5

mbuild/formats/lammpsdata.py

@@ -37,6 +37,7 @@ def write_lammpsdata(
     sigma_conversion_factor=None,
     epsilon_conversion_factor=None,
     mass_conversion_factor=None,
+    charge_conversion_factor=True,
     zero_dihedral_weighting_factor=False,
     moleculeID_offset=1,
 ):
@@ -52,51 +53,56 @@ def write_lammpsdata(
         ParmEd structure object
     filename : str
         Path of the output file
-    atom_style: str
+    atom_style: str, optional, default='full'
         Defines the style of atoms to be saved in a LAMMPS data file. The
         following atom styles are currently supported:
         'full', 'atomic', 'charge', 'molecular'
         see http://lammps.sandia.gov/doc/atom_style.html for more information
         on atom styles.
-    unit_style: str
-        Defines to unit style to be save in a LAMMPS data file.  Defaults to
-        'real' units. Current styles are supported: 'real', 'lj'
-        see https://lammps.sandia.gov/doc/99/units.html for more information
-        on unit styles
-    mins : list
-        minimum box dimension in x, y, z directions, nm
-    maxs : list
-        maximum box dimension in x, y, z directions, nm
-    pair_coeff_label : str
+    unit_style : str, optional, default='real'
+        Defines to unit style to be save in a LAMMPS data file. Current styles
+        are supported: 'real', 'lj', see lammps unit style documentation:
+        https://lammps.sandia.gov/doc/99/units.html for more information.
+    mins : list, optional, default=None
+        Minimum box dimension in x, y, z directions, nm
+    maxs : list, optional, default=None
+        Maximum box dimension in x, y, z directions, nm
+    pair_coeff_label : str, optional, default=None
         Provide a custom label to the pair_coeffs section in the lammps data
-        file. Defaults to None, which means a suitable default will be chosen.
-    detect_forcefield_style: boolean
+        file. A value of None means a suitable default will be chosen.
+    detect_forcefield_style : bool, optional, default=True
         If True, format lammpsdata parameters based on the contents of
         the parmed Structure
-    use_urey_bradleys: boolean
+    use_urey_bradleys : bool, optional, default=False
         If True, will treat angles as CHARMM-style angles with urey bradley
         terms while looking for `structure.urey_bradleys`
-    use_rb_torsions:
+    use_rb_torsions : bool, optional, default=True
         If True, will treat dihedrals OPLS-style torsions while looking for
         `structure.rb_torsions`
-    use_dihedrals:
+    use_dihedrals : bool, optional, default=False
         If True, will treat dihedrals as CHARMM-style dihedrals while looking
         for `structure.dihedrals`
-    zero_dihedral_weighting_factor:
+    zero_dihedral_weighting_factor : bool, optional, default=False
         If True, will set weighting parameter to zero in CHARMM-style dihedrals.
         This should be True if the CHARMM dihedral style is used in non-CHARMM forcefields.
-    sigma_conversion_factor: None, float
+    sigma_conversion_factor : float, optional, default=None
         If unit_style is set to 'lj', then sigma conversion factor is used to non-dimensionalize.
-        Assume to be in units of nm. Default is None. If None, will take the largest sigma value in
+        Assume to be in units of nm. If None, will take the largest sigma value in
         the structure.atoms.sigma values.
-    epsilon_conversion_factor: None, float
+    epsilon_conversion_factor : float, optional, default=None
         If unit_style is set to 'lj', then epsilon conversion factor is used to non-dimensionalize.
-        Assume to be in units of kCal/mol. Default is None. If None, will take the largest epsilon value in
+        Assume to be in units of kCal/mol. If None, will take the largest epsilon value in
         the structure.atoms.epsilon values.
-    mass_conversion_factor: None, float
+    mass_conversion_factor : float, optional, default=None
         If unit_style is set to 'lj', then mass conversion factor is used to non-dimensionalize.
-        Assume to be in units of amu. Default is None. If None, will take the largest mass value in
+        Assume to be in units of amu. If None, will take the largest mass value in
         the structure.atoms.mass values.
+    charge_conversion_factor : bool, optional, default=True
+        If unit_style is set to 'lj', then charge conversion factor may or may not be used to
+        non-dimensionalize. Assume to be in elementary charge units. If ``True``, the charges
+        are scaled by ``np.sqrt(4*np.pi()*eps_0*sigma_conversion_factor*epsilon_conversion_factor)``.
+        If ``False``, the charges are not scaled and the user must be wary to choose the dielectric
+        constant properly, which may be more convenient to implement an implicit solvent.
     moleculeID_offset : int , optional, default=1
         Since LAMMPS treats the MoleculeID as an additional set of information
         to identify what molecule an atom belongs to, this currently
@@ -287,16 +293,17 @@ def write_lammpsdata(
         )
         # Convert coordinates and charges to LJ units
         xyz = xyz / sigma_conversion_factor
-        charges = (charges * ELEM_TO_COUL) / np.sqrt(
-            4
-            * np.pi
-            * sigma_conversion_factor
-            * NM_TO_ANG**-1
-            * epsilon_conversion_factor
-            * KCAL_TO_KJ
-            * epsilon_0
-            * 10**-6
-        )
+        if charge_conversion_factor:
+            charges = (charges * ELEM_TO_COUL) / np.sqrt(
+                4
+                * np.pi
+                * sigma_conversion_factor
+                * NM_TO_ANG**-1
+                * epsilon_conversion_factor
+                * KCAL_TO_KJ
+                * epsilon_0
+                * 10**-6
+            )
         charges[np.isinf(charges)] = 0
     else:
         sigma_conversion_factor = 1