Created
March 2, 2020 16:09
-
-
Save uppittu11/c421280913edee7e3a56ce09a7d3c2d3 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import mbuild as mb\n", | |
"\n", | |
"import gmso\n", | |
"from gmso import ForceField\n", | |
"from gmso.external import convert_mbuild\n", | |
"from gmso.formats.gro import write_gro\n", | |
"from gmso.formats.top import write_top" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Write out the SPC/E water mol2 file\n", | |
"\n", | |
"with open(\"spce.mol2\", \"w\") as f:\n", | |
" f.write(\n", | |
"\"\"\"@<TRIPOS>MOLECULE\n", | |
"RES\n", | |
"3 0 1 0 1\n", | |
"SMALL\n", | |
"NO_CHARGES\n", | |
"@<TRIPOS>CRYSIN\n", | |
" 3.0130 3.0130 3.0130 90.0000 90.0000 90.0000 1 1\n", | |
"@<TRIPOS>ATOM\n", | |
" 1 O 0.0000 0.0000 0.0000 O 1 RES \n", | |
" 2 H -0.6126 -0.7357 0.0000 H 1 RES \n", | |
" 3 H -0.5469 0.7762 0.0000 H 1 RES \n", | |
"@<TRIPOS>BOND\n", | |
" 1 1 2 1\n", | |
" 2 1 3 1\n", | |
"@<TRIPOS>SUBSTRUCTURE\n", | |
" 1 RES 1 RESIDUE 0 **** ROOT 0\"\"\")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Load in single water (SPC/E) structure\n", | |
"\n", | |
"water = mb.load(\"spce.mol2\")\n", | |
"water = water.children[0]\n", | |
"water.name = \"water\"\n", | |
"\n", | |
"# element_map = which site name corresponds to which atom_type name. \n", | |
"# In the future atomtyping will be done through foyer. \n", | |
"element_map = {\"O\": \"opls_116\",\n", | |
" \"H\": \"opls_117\"}" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Fill a box with 100 water molecule\n", | |
"\n", | |
"water_box = mb.fill_box(\n", | |
" compound=water, \n", | |
" n_compounds=100,\n", | |
" density=1000,\n", | |
")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Write out the SPC/E water xml file\n", | |
"\n", | |
"with open(\"spce.xml\", \"w\") as f:\n", | |
" f.write(\n", | |
"\"\"\"<?xml version='1.0' encoding='UTF-8'?>\n", | |
"<ForceField name=\"spce\" version=\"0.0.1\">\n", | |
" <FFMetaData>\n", | |
" <Units energy=\"kcal/mol\" mass=\"amu\" charge=\"elementary_charge\" distance=\"nm\"/>\n", | |
" </FFMetaData>\n", | |
" <AtomTypes expression=\"4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)\">\n", | |
" <ParametersUnitDef parameter=\"epsilon\" unit=\"kJ/mol\"/>\n", | |
" <ParametersUnitDef parameter=\"sigma\" unit=\"nm\"/>\n", | |
" <AtomType name=\"opls_116\" atomclass=\"OW\" element=\"O\" charge=\"-0.8476\" mass=\"15.999\" definition=\"O\">\n", | |
" <Parameters>\n", | |
" <Parameter name=\"epsilon\" value=\"0.650194\"/>\n", | |
" <Parameter name=\"sigma\" value=\"0.316557\"/>\n", | |
" </Parameters>\n", | |
" </AtomType>\n", | |
" <AtomType name=\"opls_117\" atomclass=\"HW\" element=\"H\" charge=\"0.4238\" mass=\"1.00784\" definition=\"H\">\n", | |
" <Parameters>\n", | |
" <Parameter name=\"epsilon\" value=\"0.0\"/>\n", | |
" <Parameter name=\"sigma\" value=\"0.1\"/>\n", | |
" </Parameters>\n", | |
" </AtomType>\n", | |
" </AtomTypes>\n", | |
" <BondTypes expression=\"0.5 * k * (r-r_eq)**2\">\n", | |
" <ParametersUnitDef parameter=\"k\" unit=\"kJ/mol/nm**2\"/>\n", | |
" <ParametersUnitDef parameter=\"r_eq\" unit=\"nm\"/>\n", | |
" <BondType name=\"BondType-Harmonic-1\" type1=\"opls_116\" type2=\"opls_117\">\n", | |
" <Parameters>\n", | |
" <Parameter name=\"k\" value=\"345000.0\"/>\n", | |
" <Parameter name=\"r_eq\" value=\"0.1\"/>\n", | |
" </Parameters>\n", | |
" </BondType>\n", | |
" </BondTypes>\n", | |
" <AngleTypes expression=\"0.5 * k * (theta - theta_eq)**2\">\n", | |
" <ParametersUnitDef parameter=\"k\" unit=\"kJ/mol/degree**2\"/>\n", | |
" <ParametersUnitDef parameter=\"theta_eq\" unit=\"degree\"/>\n", | |
" <AngleType name=\"AngleType-Harmonic-1\" type1=\"opls_116\" type2=\"opls_117\" type3=\"opls_117\">\n", | |
" <Parameters>\n", | |
" <Parameter name=\"k\" value=\"383.0\"/>\n", | |
" <Parameter name=\"theta_eq\" value=\"1.91061193\"/>\n", | |
" </Parameters>\n", | |
" </AngleType>\n", | |
" </AngleTypes>\n", | |
"</ForceField>\"\"\")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Load in topology forcefield\n", | |
"\n", | |
"forcefield = ForceField(\"spce.xml\")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Generate a topology from the mbuild compound\n", | |
"\n", | |
"top = convert_mbuild.from_mbuild(water_box)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Assign atom types\n", | |
"for atom in top.sites:\n", | |
" atom.atom_type = forcefield.atom_types[element_map[atom.name]]\n", | |
" \n", | |
"# Assign bond types\n", | |
"for bond in top.bonds:\n", | |
" bond.bond_type = bond.connection_type = forcefield.bond_types[\"opls_116~opls_117\"]\n", | |
"\n", | |
"# Create angles with correct atom type and add to top\n", | |
"for subtop in top.subtops:\n", | |
" angle = gmso.core.angle.Angle(\n", | |
" connection_members=[site for site in subtop.sites],\n", | |
" name=\"opls_116~opls_117~opls_117\",\n", | |
" connection_type=forcefield.angle_types[\"opls_116~opls_117~opls_117\"]\n", | |
" )\n", | |
" top.add_connection(angle)\n", | |
"\n", | |
"top.update_topology()" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Write out gro and top files\n", | |
"\n", | |
"write_gro(top, \"system.gro\")\n", | |
"write_top(top, \"system.top\")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Lets look at the TOP file\n", | |
"\n", | |
"!cat system.top" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# Write out an MDP file for an NPT simulation\n", | |
"\n", | |
"with open(\"npt.mdp\", \"w\") as f:\n", | |
" f.write(\n", | |
"\"\"\"\n", | |
"constraints = h-bonds\n", | |
"integrator = md\n", | |
"nsteps = 1000000\n", | |
"dt = 0.001\n", | |
"\n", | |
"nstxtcout = 1000\n", | |
"nstenergy = 1000\n", | |
"nstlog = 1000\n", | |
"\n", | |
"cutoff-scheme = Verlet\n", | |
"ns_type = grid\n", | |
"nstlist = 10\n", | |
"rcoulomb = 0.5\n", | |
"rvdw = 0.5 \n", | |
"\n", | |
"coulombtype = PME\n", | |
"\n", | |
"gen_vel = no\n", | |
"\n", | |
"tcoupl = nose-hoover\n", | |
"tc-grps = System\n", | |
"tau_t = 1\n", | |
"ref_t = 300\n", | |
"\n", | |
"pcoupl = Parrinello-Rahman\n", | |
"pcoupltype = isotropic\n", | |
"tau-p = 10\n", | |
"ref-p = 1\n", | |
"compressibility = 1e-5\n", | |
"\"\"\"\n", | |
")" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"!gmx grompp -f npt.mdp -c system.gro -p system.top -o npt.tpr -maxwarn 2" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"!gmx mdrun -v -deffnm npt" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3", | |
"language": "python", | |
"name": "python3" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 3 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython3", | |
"version": "3.7.6" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 2 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment