Skip to content

Commit

Permalink
Merge pull request #114 from summeraz/docs
Browse files Browse the repository at this point in the history
Added usage examples markdown file
  • Loading branch information
ctk3b committed Apr 14, 2017
2 parents bf414cb + f2cd288 commit 9d8b006
Show file tree
Hide file tree
Showing 3 changed files with 112 additions and 0 deletions.
2 changes: 2 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ https://github.com/mosdef-hub/forcefield_template

### [SMARTS-based atomtyping](docs/smarts.md)

### [Usage examples](docs/usage_examples.md)

#### [![License](https://img.shields.io/badge/license-MIT-blue.svg)](http://opensource.org/licenses/MIT)

Various sub-portions of this library may be independently distributed under
Expand Down
73 changes: 73 additions & 0 deletions docs/usage_examples.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
### Usage Examples

Foyer supports atomtyping of both all-atom and coarse-grained molecular systems, and
also allows for separate force fields to be used to atom-type separate portions of a
molecular system.

#### Creating a box of ethane
Here we use mBuild to construct a box filled with ethane molecules and use Foyer to
atom-type the system, applying the OPLS force field, and save to run-able GROMACS
files.
```python
import mbuild as mb
from mbuild.examples import Ethane

ethane_box = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2])
ethane_box.save('ethane-box.gro')
ethane_box.save('ethane-box.top', forcefield_name='oplsaa')
```
----
#### Creating a box of coarse-grained ethane
Again we will use mBuild to construct a box filled with ethane molecules. However,
now we will model ethane using a united-atom description and apply the TraPPE force
field during atom-typing. Note how particle names are prefixed with an underscore so
that Foyer knows these particles are non-atomistic.
```python
import mbuild as mb

ethane_UA = mb.Compound()
ch3_1 = mb.Particle(name='_CH3', pos=[0, 0, 0])
ch3_2 = mb.Particle(name='_CH3', pos=[0.15, 0, 0])
ethane_UA.add([ch3_1, ch3_2])
ethane_UA.add_bond((ch3_1, ch3_2))

ethane_UA_box = mb.fill_box(ethane_UA, 100, box=[2, 2, 2])
ethane_UA_box.save('ethane-UA-box.gro')
ethane_UA_box.save('ethane-UA-box.top', forcefield_name='trappe-ua')
```
----
#### Combining force fields
In some instances, the use of multiple force fields may be desired to describe a
molecular system. For example, the user may want to use one force field for a
surface and another for a fluid in the same system. Foyer supports this
functionality by allowing the user to separately atom-type parts of a system. In
this example, we take a system featuring bulk united atom ethane above a silica
surface and apply the OPLS force field to the surface and the TraPPE force field to
the ethane. The two atomtyped Parmed structures are then combined using a simple
'\+' operator and can be saved to Gromacs files.
```python
from foyer import Forcefield
from foyer.test.utils import get_fn
import mbuild as mb
from mbuild.examples import Ethane
from mbuild.lib.atoms import H
from mbuild.lib.bulk_materials import AmorphousSilica

interface = mb.SilicaInterface(bulk_silica=AmorphousSilica())
interface = mb.Monolayer(surface=interface, chains=H(), guest_port_name='up')

box = mb.Box(mins=[0, 0, max(interface.xyz[:,2])],
maxs=interface.periodicity + [0, 0, 4])

ethane_box = mb.fill_box(compound=Ethane(), n_compounds=200, box=box)

opls = Forcefield(name='oplsaa')
opls_silica = Forcefield(forcefield_files=get_fn('opls-silica.xml'))
ethane_box = opls.apply(ethane_box)
interface = opls_silica.apply(interface)

system = interface + ethane_box

system.save('ethane-silica.gro')
system.save('ethane-silica.top')
```
37 changes: 37 additions & 0 deletions foyer/tests/files/opls-silica.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
<ForceField>
<AtomTypes>
<Type name="opls_1001" class="OS" element="O" mass="15.9994" def="[O][Si;%opls_1002]"/>
<Type name="opls_1002" class="SI" element="Si" mass="28.0855" def="[Si]"/>
<Type name="opls_1003" class="OS" element="O" mass="15.9994" def="[O;X2]([Si;%opls_1002])H" overrides="opls_1001"/>
<Type name="opls_1004" class="HO" element="H" mass="1.00800" def="[H;X1](O[Si;%opls_1002])"/>
</AtomTypes>
<HarmonicBondForce>
<Bond class1="SI" class2="OS" length="0.163" k="251040.0"/>
<Bond class1="SI" class2="OH" length="0.163" k="251040.0"/>
<Bond class1="SI" class2="CT" length="0.185" k="167360.0"/>
</HarmonicBondForce>
<HarmonicAngleForce>
<Angle class1="SI" class2="OS" class3="SI" angle="2.7053" k="397.48"/>
<Angle class1="SI" class2="OH" class3="SI" angle="2.7053" k="397.48"/>
<Angle class1="OS" class2="SI" class3="OS" angle="1.9111" k="397.48"/>
<Angle class1="OH" class2="SI" class3="OS" angle="1.9111" k="397.48"/>
<Angle class1="SI" class2="OH" class3="HO" angle="1.9111" k="397.48"/>
<Angle class1="SI" class2="CT" class3="CT" angle="2.0944" k="254.97296"/>
<Angle class1="SI" class2="CT" class3="HC" angle="1.9111" k="418.48"/>
<Angle class1="OH" class2="SI" class3="CT" angle="1.7453" k="502.18"/>
<Angle class1="OH" class2="SI" class3="OH" angle="1.9199" k="502.18"/>
<Angle class1="OS" class2="SI" class3="CT" angle="1.7453" k="502.18"/>
</HarmonicAngleForce>
<RBTorsionForce>
<Proper class1="" class2="SI" class3="OS" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
<Proper class1="" class2="SI" class3="OH" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
<Proper class1="" class2="SI" class3="CT" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
<Proper class1="SI" class2="CT" class3="CT" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
</RBTorsionForce>
<NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
<Atom type="opls_1001" charge="-0.43" sigma="0.3" epsilon="0.71128"/>
<Atom type="opls_1002" charge="0.86" sigma="0.4" epsilon="0.4184"/>
<Atom type="opls_1003" charge="-0.215" sigma="0.3" epsilon="0.71128"/>
<Atom type="opls_1004" charge="0.215" sigma="0.0" epsilon="0.0"/>
</NonbondedForce>
</ForceField>

0 comments on commit 9d8b006

Please sign in to comment.