Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Map forces onto beads #80

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open

Map forces onto beads #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
5766156
adding comments for things to do
StephMcCallum Jun 22, 2024
00a4681
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 22, 2024
365d643
Adding force list
StephMcCallum Jun 22, 2024
1a63d1e
Adding force list
StephMcCallum Jun 22, 2024
8b51bdd
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 22, 2024
83c887a
Adding test velocity calculation
StephMcCallum Jun 25, 2024
eb743e7
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 26, 2024
b0f80dc
Adding force calculation and test file. Need to update method of addi…
StephMcCallum Aug 8, 2024
2a4194f
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Aug 8, 2024
6d35f1a
Changing axis of reduce function and adding squeeze function for a te…
StephMcCallum Aug 12, 2024
03a04b4
changing axis for reduce function and adding squeeze
StephMcCallum Aug 12, 2024
660bf1c
Updating net force calculation
StephMcCallum Aug 13, 2024
13eda05
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Aug 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
273 changes: 273 additions & 0 deletions Test run.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,273 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "86a0a551-a4df-4664-b230-63c0b6577256",
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/stephaniemccallum/miniforge3/envs/grits/lib/python3.12/site-packages/mdtraj/formats/__init__.py:13: DeprecationWarning: 'xdrlib' is deprecated and slated for removal in Python 3.13\n",
" from mdtraj.formats.trr import TRRTrajectoryFile\n",
"/Users/stephaniemccallum/miniforge3/envs/grits/lib/python3.12/site-packages/foyer/forcefield.py:33: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html\n",
" from pkg_resources import iter_entry_points, resource_filename\n",
"/Users/stephaniemccallum/miniforge3/envs/grits/lib/python3.12/site-packages/pkg_resources/__init__.py:3144: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('google')`.\n",
"Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages\n",
" declare_namespace(pkg)\n",
"/Users/stephaniemccallum/miniforge3/envs/grits/lib/python3.12/site-packages/lark/utils.py:163: DeprecationWarning: module 'sre_parse' is deprecated\n",
" import sre_parse\n",
"/Users/stephaniemccallum/miniforge3/envs/grits/lib/python3.12/site-packages/lark/utils.py:164: DeprecationWarning: module 'sre_constants' is deprecated\n",
" import sre_constants\n",
"/Users/stephaniemccallum/miniforge3/envs/grits/lib/python3.12/site-packages/mbuild/packing.py:23: DeprecationWarning: Use shutil.which instead of find_executable\n",
" PACKMOL = find_executable(\"packmol\")\n"
]
}
],
"source": [
"import numpy as np\n",
"import mbuild as mb\n",
"\n",
"from grits import CG_System\n",
"from grits.utils import amber_dict"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "52fd0a7c-35af-43b2-874f-a6b52fc03d65",
"metadata": {},
"outputs": [],
"source": [
"PPS_smiles = \"c1ccc(S)cc1\"\n",
"a = mb.load(PPS_smiles, smiles=True)\n",
"#a.visualize().show()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "d6285f42-272c-4206-b738-fa3cb3688626",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Added 6 hydrogens.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/stephaniemccallum/Desktop/cmelab/forks/grits/grits/coarsegrain.py:192: UserWarning: Some atoms have been left out of coarse-graining!\n",
" warn(f\"Some atoms have been left out of coarse-graining!\")\n"
]
}
],
"source": [
"gsdfile = \"/Users/stephaniemccallum/Desktop/cmelab/forks/grits/pps_trajectory_kT_2.0.gsd\"\n",
"system = CG_System(\n",
" gsdfile,\n",
" add_hydrogens = True,\n",
" beads={\"_A\": PPS_smiles},\n",
" conversion_dict=amber_dict\n",
")\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "9660225b-93e6-49b8-9d31-5786ae97e187",
"metadata": {},
"outputs": [],
"source": [
"cg_gsd = \"PPS-cg.gsd\"\n",
"system.save(cg_gsd)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "90ecc7df-c73d-489d-8053-143459f91313",
"metadata": {},
"outputs": [],
"source": [
"import gsd.hoomd\n",
"\n",
"traj = gsd.hoomd.open(gsdfile)\n",
"snap = traj[-1]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "7a434d67-2d79-4ebf-b400-709a5baa8c66",
"metadata": {},
"outputs": [],
"source": [
"inds = [0,1,2,3]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "5474c6f6-3fad-49b8-90cd-6142b870b379",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-0.58440879, -0.72004115, -0.20885618])"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"\n",
"np.add.reduce(snap.log[\"particles/md/pair/LJ/forces\"][inds])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "4d0906fe-0023-4b98-9cdb-97f133843a8b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-0.59678146 -0.7094508 -0.20883539]\n",
" [ 0. 0. 0. ]\n",
" [ 0. 0. 0. ]\n",
" ...\n",
" [-0.23439395 -0.89780405 -0.28628996]\n",
" [ 0. 0. 0. ]\n",
" [ 0. 0. 0. ]]\n"
]
}
],
"source": [
"print(snap.log[\"particles/md/pair/LJ/forces\"])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "36fbdb7b-7449-4cbe-bc1e-40f5f158b1dd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.46751893, -0.4557764 , 1.5517997 ],\n",
" [ 1.2661922 , 2.2942452 , 2.3974137 ],\n",
" [ 0.99841857, -0.54276377, -0.38419962],\n",
" ...,\n",
" [-0.16005144, 0.94953936, -1.093895 ],\n",
" [ 0.6594045 , -1.5998036 , -1.1851496 ],\n",
" [-0.8315832 , 1.2608055 , 0.19617926]], dtype=float32)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"snap.particles.velocity"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "6ddec25f-8367-410a-91f1-e18082affb6f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['__class__',\n",
" '__delattr__',\n",
" '__dict__',\n",
" '__dir__',\n",
" '__doc__',\n",
" '__eq__',\n",
" '__format__',\n",
" '__ge__',\n",
" '__getattribute__',\n",
" '__getstate__',\n",
" '__gt__',\n",
" '__hash__',\n",
" '__init__',\n",
" '__init_subclass__',\n",
" '__le__',\n",
" '__lt__',\n",
" '__module__',\n",
" '__ne__',\n",
" '__new__',\n",
" '__reduce__',\n",
" '__reduce_ex__',\n",
" '__repr__',\n",
" '__setattr__',\n",
" '__sizeof__',\n",
" '__slotnames__',\n",
" '__str__',\n",
" '__subclasshook__',\n",
" '__weakref__',\n",
" '_valid_state',\n",
" 'angles',\n",
" 'bonds',\n",
" 'configuration',\n",
" 'constraints',\n",
" 'dihedrals',\n",
" 'impropers',\n",
" 'log',\n",
" 'pairs',\n",
" 'particles',\n",
" 'state',\n",
" 'validate']"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dir(snap)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
52 changes: 51 additions & 1 deletion grits/coarsegrain.py
View file
Open in desktop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the additons so far look good. Now, you'll have to add the coarse-grained velocities to the new gsd frame that GRiTS is making.

You can see what is done for other properties around lines 780:

                new_snap.particles.N = len(typeid)
                new_snap.particles.position = position
                new_snap.particles.image = images
                new_snap.particles.typeid = typeid.astype(int)
                new_snap.particles.types = types
                new_snap.particles.mass = mass

So, you can define the new_snap.particles.velocity property with the velocity list being populated.

Original file line number Diff line number Diff line change
Expand Up @@ -727,6 +727,14 @@ def save(self, cg_gsdfile, start=0, stop=None, stride=1):
new_snap = gsd.hoomd.Frame()
position = []
mass = []
# make an empty lists for forces here
traj_lj_forces = []
traj_bond_forces = []
traj_angle_forces = []
traj_dihedral_forces = []
net_force = []
# make an empty list for velocity
velocity = []
orientation = [] if self.aniso_beads else None
f_box = freud.Box.from_box(s.configuration.box)
unwrap_pos = f_box.unwrap(
Expand All @@ -738,6 +746,40 @@ def save(self, cg_gsdfile, start=0, stop=None, stride=1):
mass += [
sum(s.particles.mass[x]) * self.mass_scale for x in inds
]

# do the force calculation here
traj_lj_forces.append(
np.add.reduce(
s.log["particles/md/pair/LJ/forces"][inds], 1
)
)
traj_bond_forces.append(
np.add.reduce(
s.log["particles/md/pair/LJ/forces"][inds], 1
)
)
traj_angle_forces.append(
np.add.reduce(
s.log["particles/md/pair/LJ/forces"][inds], 1
)
)
traj_dihedral_forces.append(
np.add.reduce(
s.log["particles/md/pair/LJ/forces"][inds], 1
)
)
traj_lj_forces = np.squeeze(traj_lj_forces, axis=0)
traj_bond_forces = np.squeeze(traj_bond_forces, axis=0)
traj_angle_forces = np.squeeze(traj_angle_forces, axis=0)
traj_dihedral_forces = np.squeeze(
traj_dihedral_forces, axis=0
)

# do the velocity calculation here
velocity += [
np.mean(s.particles.velocity[x], axis=0) for x in inds
]

if self.aniso_beads:
for x in inds:
masses = s.particles.mass[x] * self.mass_scale
Expand All @@ -757,7 +799,13 @@ def save(self, cg_gsdfile, start=0, stop=None, stride=1):
position = np.vstack(position)
images = f_box.get_images(position)
position = f_box.wrap(position)

arr = [
np.array(traj_lj_forces),
np.array(traj_bond_forces),
np.array(traj_angle_forces),
np.array(traj_dihedral_forces),
]
net_force = np.add.reduce(arr)
new_snap.configuration.box = s.configuration.box
new_snap.configuration.step = s.configuration.step
new_snap.particles.N = len(typeid)
Expand All @@ -766,6 +814,8 @@ def save(self, cg_gsdfile, start=0, stop=None, stride=1):
new_snap.particles.typeid = typeid.astype(int)
new_snap.particles.types = types
new_snap.particles.mass = mass
new_snap.particles.velocity = velocity
new_snap.log["net_force"] = np.array(net_force)

if N_bonds > 0:
new_snap.bonds.N = N_bonds
Expand Down
Loading