diff --git a/Test run.ipynb b/Test run.ipynb new file mode 100644 index 0000000..d5bfbff --- /dev/null +++ b/Test run.ipynb @@ -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 +} diff --git a/grits/coarsegrain.py b/grits/coarsegrain.py index 7d7fd2f..1444c07 100644 --- a/grits/coarsegrain.py +++ b/grits/coarsegrain.py @@ -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( @@ -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 @@ -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) @@ -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