Skip to content

Commit

Permalink
Update repo
Browse files Browse the repository at this point in the history
  • Loading branch information
@gitesei
gitesei committed Dec 14, 2023
1 parent 82584cf commit 4a1a426
Show file tree
Hide file tree
Showing 72 changed files with 304,297 additions and 23,417 deletions.
234 changes: 113 additions & 121 deletions analyses.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,18 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 41,
"id": "29e13279",
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/Users/vbr805/Dropbox/lnps_manuscript/SI-LNPs-2\n"
]
}
],
"source": [
"import numpy as np\n",
"import mdtraj as md\n",
Expand Down Expand Up @@ -573,83 +581,84 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 42,
"id": "5709d5d0",
"metadata": {},
"outputs": [],
"source": [
"models_snapshots = [3.09,2.54,1.86,1.86][:]\n",
"vmd_path = 'csh /Applications/VMD_1.9.4a57-arm64-Rev12.app/Contents/MacOS/startup.command.csh'\n",
"if os.path.isfile(vmd_path):\n",
" models_snapshots = [3.09,2.54,1.86,1.86][:]\n",
"\n",
"pH_values = [3.0, 4.0, 5.0, 5.8, 6.6, 7.4, 8.0][:3]\n",
" pH_values = [3.0, 4.0, 5.0, 5.8, 6.6, 7.4, 8.0][:3]\n",
"\n",
"def writePQRAtom(line,ofh):\n",
" spl = line.split()\n",
" aname=spl[2]\n",
" if not aname[0].isdigit():\n",
" #print \"change aname\\\"%s\\\"\"%aname\n",
" aname = \" \"+aname\n",
" #print \"to \\\"%s\\\"\"%aname\n",
" if len(spl)==11:\n",
" outline=\"ATOM %5d %-4s%1s%3s %1s%4d%1s %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s%s\\n\" % (\n",
" int(spl[1]),aname,\" \",spl[3],spl[4],int(spl[5]),\" \",float(spl[6]),float(spl[7]),\n",
" float(spl[8]),float(spl[9]),float(spl[10]),\" \",\" \",\" \",\"\")\n",
" elif len(spl)==10:\n",
" outline=\"ATOM %5d %-4s%1s%3s %1s%4d%1s %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s%s\\n\" % (\n",
" int(spl[1]),aname,\" \",spl[3],\" \",int(spl[4]),\" \",float(spl[5]),float(spl[6]),\n",
" float(spl[7]),float(spl[8]),float(spl[9]),\" \",\" \",\" \",\"\")\n",
" else:\n",
" raise ValueError(\"The format of this pqr is not as expected\")\n",
" ofh.write(outline)\n",
" \n",
"def hex_to_vmd(color):\n",
" h = mpl.colors.to_hex(color).lstrip('#')\n",
" rgb = np.array([ int(h[i:i+2], 16) for i in (0, 2 ,4)])\n",
" return ' '.join(['{:1.2f}'.format(i/255.) for i in rgb])\n",
" def writePQRAtom(line,ofh):\n",
" spl = line.split()\n",
" aname=spl[2]\n",
" if not aname[0].isdigit():\n",
" #print \"change aname\\\"%s\\\"\"%aname\n",
" aname = \" \"+aname\n",
" #print \"to \\\"%s\\\"\"%aname\n",
" if len(spl)==11:\n",
" outline=\"ATOM %5d %-4s%1s%3s %1s%4d%1s %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s%s\\n\" % (\n",
" int(spl[1]),aname,\" \",spl[3],spl[4],int(spl[5]),\" \",float(spl[6]),float(spl[7]),\n",
" float(spl[8]),float(spl[9]),float(spl[10]),\" \",\" \",\" \",\"\")\n",
" elif len(spl)==10:\n",
" outline=\"ATOM %5d %-4s%1s%3s %1s%4d%1s %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s%s\\n\" % (\n",
" int(spl[1]),aname,\" \",spl[3],\" \",int(spl[4]),\" \",float(spl[5]),float(spl[6]),\n",
" float(spl[7]),float(spl[8]),float(spl[9]),\" \",\" \",\" \",\"\")\n",
" else:\n",
" raise ValueError(\"The format of this pqr is not as expected\")\n",
" ofh.write(outline)\n",
"\n",
"for folder in ['lnp_w_rna_310','lnp_empty_310']:\n",
" path = f'data/{folder:s}'+'/CPP{:.2f}_{:.2f}_{:.1f}'\n",
" for model,pKa in zip(models_snapshots,[9.41,9.41,9.41,8.46]):\n",
" l = []\n",
" for pH in pH_values:\n",
" current_path = path.format(model,pKa,pH)\n",
" ifname = current_path+'/cuboid_1000000.pqr'\n",
" ofname = current_path+'/cuboid_1000000.pdb'\n",
" with open(ifname,\"r\") as ifhandle:\n",
" with open(ofname,\"w\") as ofhandle:\n",
" for line in ifhandle:\n",
" if line.startswith(\"ATOM\"):\n",
" writePQRAtom(line,ofhandle)\n",
" else:\n",
" ofhandle.write(line)\n",
" s = md.load_pdb(ofname)\n",
" index_1 = s.top.select('name T1')\n",
" index_2 = s.top.select('name T2')\n",
" d = md.compute_distances(s,np.c_[index_1,index_2])\n",
" l.append(d[d<3].mean())\n",
" def hex_to_vmd(color):\n",
" h = mpl.colors.to_hex(color).lstrip('#')\n",
" rgb = np.array([ int(h[i:i+2], 16) for i in (0, 2 ,4)])\n",
" return ' '.join(['{:1.2f}'.format(i/255.) for i in rgb])\n",
"\n",
"colors = ['tab:green','tab:blue','tab:orange','tab:red']\n",
"vmd_path = 'csh /Applications/VMD_1.9.4a57-arm64-Rev12.app/Contents/MacOS/startup.command.csh'\n",
" for folder in ['lnp_w_rna_310','lnp_empty_310']:\n",
" path = f'data/{folder:s}'+'/CPP{:.2f}_{:.2f}_{:.1f}'\n",
" for model,pKa in zip(models_snapshots,[9.41,9.41,9.41,8.46]):\n",
" l = []\n",
" for pH in pH_values:\n",
" current_path = path.format(model,pKa,pH)\n",
" ifname = current_path+'/cuboid_1000000.pqr'\n",
" ofname = current_path+'/cuboid_1000000.pdb'\n",
" with open(ifname,\"r\") as ifhandle:\n",
" with open(ofname,\"w\") as ofhandle:\n",
" for line in ifhandle:\n",
" if line.startswith(\"ATOM\"):\n",
" writePQRAtom(line,ofhandle)\n",
" else:\n",
" ofhandle.write(line)\n",
" s = md.load_pdb(ofname)\n",
" index_1 = s.top.select('name T1')\n",
" index_2 = s.top.select('name T2')\n",
" d = md.compute_distances(s,np.c_[index_1,index_2])\n",
" l.append(d[d<3].mean())\n",
"\n",
" colors = ['tab:green','tab:blue','tab:orange','tab:red']\n",
"\n",
"for folder in ['lnp_w_rna_310','lnp_empty_310']:\n",
" path = f'data/{folder:s}'+'/CPP{:.2f}_{:.2f}_{:.1f}'\n",
" script_path = f'../snapshot.vmd'\n",
" for pH in pH_values:\n",
" for model,pKa,color in zip(models_snapshots,[9.41,9.41,9.41,8.46],colors):\n",
" current_path = path.format(model,pKa,pH)\n",
" %cd {current_path}\n",
" # Read in the file\n",
" with open(script_path, 'r') as file :\n",
" filedata = file.read()\n",
" filedata = filedata.replace('new_box_color', hex_to_vmd(color))\n",
" if folder == 'lnp_empty_310':\n",
" half_box_length_y = md.load('confout.gro').unitcell_lengths[0][0]*np.sqrt(3)*10/2\n",
" filedata = filedata.replace('half_box_length_y', f'{half_box_length_y:f}')\n",
" with open('snapshot.vmd', 'w') as file:\n",
" file.write(filedata)\n",
" print(f'../{current_path.split(\"/\")[-1]:s}.tga')\n",
" !{vmd_path} -e snapshot.vmd -dispdev text\n",
" !cp snapshot.tga ../{current_path.split('/')[-1]}.tga\n",
" %cd -"
" for folder in ['lnp_w_rna_310','lnp_empty_310']:\n",
" path = f'data/{folder:s}'+'/CPP{:.2f}_{:.2f}_{:.1f}'\n",
" script_path = f'../snapshot.vmd'\n",
" for pH in pH_values:\n",
" for model,pKa,color in zip(models_snapshots,[9.41,9.41,9.41,8.46],colors):\n",
" current_path = path.format(model,pKa,pH)\n",
" %cd {current_path}\n",
" # Read in the file\n",
" with open(script_path, 'r') as file :\n",
" filedata = file.read()\n",
" filedata = filedata.replace('new_box_color', hex_to_vmd(color))\n",
" if folder == 'lnp_empty_310':\n",
" half_box_length_y = md.load('confout.gro').unitcell_lengths[0][0]*np.sqrt(3)*10/2\n",
" filedata = filedata.replace('half_box_length_y', f'{half_box_length_y:f}')\n",
" with open('snapshot.vmd', 'w') as file:\n",
" file.write(filedata)\n",
" print(f'../{current_path.split(\"/\")[-1]:s}.tga')\n",
" !{vmd_path} -e snapshot.vmd -dispdev text\n",
" !cp snapshot.tga ../{current_path.split('/')[-1]}.tga\n",
" %cd -"
]
},
{
Expand Down Expand Up @@ -1388,64 +1397,47 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 43,
"id": "16093987",
"metadata": {},
"outputs": [],
"source": [
"pH_values = [3.0]\n",
"\n",
"vmd_path = 'csh /Applications/VMD_1.9.4a57-arm64-Rev12.app/Contents/MacOS/startup.command.csh'\n",
"if os.path.isfile(vmd_path):\n",
"\n",
"Rcyls_298 = {}\n",
" pH_values = [3.0]\n",
"\n",
"for folder in ['lnp_empty_298']: \n",
" path = f'data/{folder:s}'+'/CPP{:.2f}_{:.2f}_{:.1f}'\n",
" script_path = f'../snapshot.vmd'\n",
" for pH in pH_values:\n",
" for model,pKa,y_sel,out_sel in zip([1.86,3.09],[8.46,9.41],\n",
" ['y<40','(y>40 and y<80)'],\n",
" ['y>-25 and y<40','y>30 and y<85']):\n",
" Rcyls_298[model] = []\n",
" print('MODEL',model)\n",
" for frame_number in range(1000,2001000,1000):\n",
" current_path = path.format(model,pKa,pH)\n",
" %cd {current_path}\n",
" # Read in the file\n",
" with open(script_path, 'r') as file :\n",
" filedata = file.read()\n",
" filedata = filedata.replace('frame_number', f'{frame_number:d}')\n",
" filedata = filedata.replace('y_sel', f'and {y_sel:s}')\n",
" filedata = filedata.replace('out_sel', f'{out_sel:s}')\n",
" with open('snapshot.vmd', 'w') as file:\n",
" file.write(filedata)\n",
" !{vmd_path} -e snapshot.vmd -dispdev text > /dev/null 2>&1\n",
" s = md.load(f'{frame_number:d}.pdb')\n",
" s_hd = s.atom_slice(s.top.select('name HD or name HHD'))\n",
" center = md.compute_center_of_geometry(s_hd)[0]\n",
" Rcyls_298[model].append(np.sqrt((s_hd.xyz[0,:,0]-center[0])**2+(s_hd.xyz[0,:,1]-center[1])**2).mean()*10)\n",
" %cd -"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c7eabfd3",
"metadata": {},
"outputs": [],
"source": [
"np.savetxt(f'data/lnp_empty_298/CPP1.86_8.46_3.0/Rcyl.dat',\n",
" np.c_[np.arange(len(Rcyls_298[1.86]))*1000,Rcyls_298[1.86]])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "508a300e",
"metadata": {},
"outputs": [],
"source": [
"np.savetxt(f'data/lnp_empty_298/CPP3.09_9.41_3.0/Rcyl.dat',\n",
" Rcyls_298 = {}\n",
"\n",
" for folder in ['lnp_empty_298']: \n",
" path = f'data/{folder:s}'+'/CPP{:.2f}_{:.2f}_{:.1f}'\n",
" script_path = f'../snapshot.vmd'\n",
" for pH in pH_values:\n",
" for model,pKa,y_sel,out_sel in zip([1.86,3.09],[8.46,9.41],\n",
" ['y<40','(y>40 and y<80)'],\n",
" ['y>-25 and y<40','y>30 and y<85']):\n",
" Rcyls_298[model] = []\n",
" print('MODEL',model)\n",
" for frame_number in range(1000,2001000,1000):\n",
" current_path = path.format(model,pKa,pH)\n",
" %cd {current_path}\n",
" # Read in the file\n",
" with open(script_path, 'r') as file :\n",
" filedata = file.read()\n",
" filedata = filedata.replace('frame_number', f'{frame_number:d}')\n",
" filedata = filedata.replace('y_sel', f'and {y_sel:s}')\n",
" filedata = filedata.replace('out_sel', f'{out_sel:s}')\n",
" with open('snapshot.vmd', 'w') as file:\n",
" file.write(filedata)\n",
" !{vmd_path} -e snapshot.vmd -dispdev text > /dev/null 2>&1\n",
" s = md.load(f'{frame_number:d}.pdb')\n",
" s_hd = s.atom_slice(s.top.select('name HD or name HHD'))\n",
" center = md.compute_center_of_geometry(s_hd)[0]\n",
" Rcyls_298[model].append(np.sqrt((s_hd.xyz[0,:,0]-center[0])**2+(s_hd.xyz[0,:,1]-center[1])**2).mean()*10)\n",
" %cd -\n",
" np.savetxt(f'data/lnp_empty_298/CPP1.86_8.46_3.0/Rcyl.dat',\n",
" np.c_[np.arange(len(Rcyls_298[1.86]))*1000,Rcyls_298[1.86]])\n",
" np.savetxt(f'data/lnp_empty_298/CPP3.09_9.41_3.0/Rcyl.dat',\n",
" np.c_[np.arange(len(Rcyls_298[3.09]))*1000,Rcyls_298[3.09]])"
]
}
Expand Down
Loading

0 comments on commit 4a1a426

Please sign in to comment.