Acellera / htmd

HTMD: Programming Environment for Molecular Discovery
https://software.acellera.com/docs/latest/htmd/index.html
Other
253 stars 58 forks source link

Building a membrane with more than three different lipids #1045

Open Saleh-OM4R opened 1 year ago

Saleh-OM4R commented 1 year ago

Hi,

I'm trying to build a membrane with more than three different lipids. I managed to add the new lipids to lipiddb.csv and add the pdb files in the share folder, I also could build a membrane with the new lipids as long as I don't exceed three lipids. Please see below the code I used and the error that I got.

HTMD Version (2.2.0)

Input


##Import
from htmd.membranebuilder.build_membrane import listLipids, buildMembrane

# List all available lipids
listLipids()
# Define the dimensions of the membrane in the x and y axis in units of Angstrom
dimensions = [100, 100]
# Define the upper and lower layer lipid composition and ratio
#ratioupper = {'popc':20,'pope':3,'popa':10}
#ratiolower = {'popc':7,'pope':14,'chl1':29}
ratioupper = {'popa':1,'nsm':2,'ssm':2,'chl1':2}
ratiolower = {'popa':1,'nsm':2,'ssm':2,'chl1':2}
# Build the membrane
memb = buildMembrane(dimensions, ratioupper, ratiolower)
# Visualize the membrane
#memb.view()

Error

---- Lipids list: /home/saleh/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/share/membranebuilder/lipids/ ----
-  chl1
-  nsm
-  pape
-  paps
-  plpc
-  popa
-  popc
-  pope
-  popi
-  pops
-  psm
-  ssm
* Lipid DB file: /home/saleh/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/share/membranebuilder/lipids/lipiddb.csv
/home/saleh/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/membranebuilder/build_membrane.py:164: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  cycles = np.array(nx.cycle_basis(G))
/home/saleh/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/membranebuilder/build_membrane.py:164: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  cycles = np.array(nx.cycle_basis(G))
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In [12], line 14
     12 ratiolower = {'popa':1,'nsm':2,'ssm':2,'chl1':2}
     13 # Build the membrane
---> 14 memb = buildMembrane(dimensions, ratioupper, ratiolower)

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/membranebuilder/build_membrane.py:284, in buildMembrane(xysize, ratioupper, ratiolower, waterbuff, minimplatform, equilibrate, equilplatform, outdir, lipidf)
    281 lipids = _createLipids(ratioupper, area, lipiddb, files, leaflet="upper")
    282 lipids += _createLipids(ratiolower, area, lipiddb, files, leaflet="lower")
--> 284 _setPositionsLJSim(xysize, [ll for ll in lipids if ll.xyz[2] > 0])
    285 _setPositionsLJSim(xysize, [ll for ll in lipids if ll.xyz[2] < 0])
    287 _findNeighbours(lipids, xysize)

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/membranebuilder/build_membrane.py:128, in _setPositionsLJSim(width, lipids)
    125 resnames = [ll.resname for ll in lipids]
    127 cutoff = min(np.min(width) / 2, 3 * np.max(sigmas))
--> 128 pos = distributeLipids(width + [2 * cutoff], resnames, sigmas, cutoff=cutoff)
    129 for i in range(len(lipids)):
    130     lipids[i].xyz[:2] = pos[i, :2]

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/htmd/membranebuilder/ljfluid.py:136, in distributeLipids(boxsize, resnames, sigmas, cutoff, mass, epsilon, switch_width)
    134 _, idx = np.unique(resnames, return_inverse=True)
    135 for i in idx:
--> 136     element = app.Element.getBySymbol(elems[i])
    137     residue = topology.addResidue(elems[i], chain)
    138     topology.addAtom(elems[i], element, residue)

IndexError: list index out of range

Regards,

stefdoerr commented 1 year ago

Oof nice catch. I wrote this code like 6 years ago as an experiment and there are apparently some hacks there which limit it to 3 lipids. I'll see what I can do about it.

stefdoerr commented 1 year ago

Ok I believe I have fixed it now. Would you mind testing the latest master branch to see if it works for you? Otherwise you can wait for the next HTMD release but it might take a bit.