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.

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

RuntimeError (adaptive MD for RNA): when using a dihedral with an atom name containing a single quote #1056

Open
vas2201 opened this issue Mar 21, 2023 · 13 comments

Comments

@vas2201
Copy link

vas2201 commented Mar 21, 2023

I encountered a RuntimeError when running the ad.run() function using a dihedral with an atom name containing a single quote (e.g., "O3'"). Here's the code I used:

atom1_alpha_Dihang_residue_6 = {'name': "O3'", 'resid': 5, 'segid': 'RNAA', 'chain': ' '}
atom2_alpha_Dihang_residue_6 = {'name': 'P', 'resid': 6, 'segid': 'RNAA', 'chain': ' '}
atom3_alpha_Dihang_residue_6 = {'name': 'O5', 'resid': 6, 'segid': 'RNAA', 'chain': ' '}
atom4_alpha_Dihang_residue_6 = {'name': 'C5', 'resid': 6, 'segid': 'RNAA', 'chain': ' '}
dih_alpha_6 = [Dihedral(atom1_alpha_Dihang_residue_6, atom2_alpha_Dihang_residue_6, atom3_alpha_Dihang_residue_6, atom4_alpha_Dihang_residue_6)]

all_dih = dih_alpha_6
print(all_dih)

# The output is as follows:
# name	resid	insertion	chain	segid
# O3'	5			 	RNAA
# P	6			 	RNAA
# O5	6			 	RNAA
# C5	6			 	RNAA

**************************************************************************
When I proceeded to run ad.run(), I got the following error:
RuntimeError: Expected one atom from atomselection {'name': "O3'", 'resid': 5, 'segid': 'RNAA', 'chain': ' ', 'insertion': ''}. Got 0 instead.

I believe the issue might be related to the atom name containing a single quote.

Please advise on this matter.
note : I uploaded notebook for quick look.
https://drive.google.com/file/d/148kbylnU8BTyF7sTz4xXRNbbQnNBVWrs/view?usp=share_link
Thank you

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[88], line 1
----> 1 ad.run()

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/adaptive/adaptive.py:196, in AdaptiveBase.run(self)
    194 # If currently running simulations are lower than nmin start new ones to reach nmax number of sims
    195 if self._running <= self.nmin and epoch < self.nepochs:
--> 196     flag = self._algorithm()
    197     if flag is False:
    198         self._unsetLock()

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/adaptive/adaptiverun.py:208, in AdaptiveMD._algorithm(self)
    207 def _algorithm(self):
--> 208     data = self._getData(self._getSimlist())
    209     if not self._checkNFrames(data):
    210         return False

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/adaptive/adaptiverun.py:246, in AdaptiveMD._getData(self, sims)
    241 # if self.contactsym is not None:
    242 #    contactSymmetry(data, self.contactsym)
    244 if self.ticadim > 0:
    245     # tica = TICA(metr, int(max(2, np.ceil(self.ticalag))))  # gianni: without project it was tooooo slow
--> 246     data = metr.project()
    247     data.dropTraj()  # Drop before TICA to avoid broken trajectories
    248     # 1 < ticalag < (trajLen / 2)

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/projections/metric.py:159, in Metric.project(self, njobs)
    157     for proj in self.projectionlist:
    158         if isinstance(proj, Projection):
--> 159             proj._setCache(uqMol)
    160 else:
    161     logger.warning(
    162         "Cannot calculate description of dimensions due to different topology files for each trajectory."
    163     )

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/moleculekit/projections/projection.py:31, in Projection._setCache(self, mol)
     30 def _setCache(self, mol):
---> 31     resdict = self._calculateMolProp(mol)
     32     self._cache.update(resdict)

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/moleculekit/projections/metricdihedral.py:839, in MetricDihedral._calculateMolProp(self, mol, props)
    836 else:
    837     dihedrals = ensurelist(self._dihedrals)
--> 839 res["dihedrals"] = Dihedral.dihedralsToIndexes(mol, dihedrals, protsel)
    840 return res

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/moleculekit/projections/metricdihedral.py:136, in Dihedral.dihedralsToIndexes(mol, dihedrals, sel)
    134     atomsel = atomsel & selatoms
    135     if np.sum(atomsel) != 1:
--> 136         raise RuntimeError(
    137             "Expected one atom from atomselection {}. Got {} instead.".format(
    138                 a, np.sum(atomsel)
    139             )
    140         )
    141     idx.append(np.where(atomsel)[0][0])
    142 indexes.append(idx)

RuntimeError: Expected one atom from atomselection {'name': "O3'", 'resid': 5, 'segid': 'RNAA', 'chain': ' ', 'insertion': ''}. Got 0 instead.
@vas2201
Copy link
Author

vas2201 commented Mar 22, 2023

FYI ...
If I use single quote around O3' atom as follows ..
atom1_alpha_Dihang_residue_6 = {'name': 'O3'', 'resid': 5, 'segid': 'RNAA', 'chain': ' '}

it is causing syntax error , that pointing out at 'resid'
Cell In[30], line 1
atom1_alpha_Dihang_residue_6 = {'name': 'O3'', 'resid': 5, 'segid': 'RNAA', 'chain': ' '}

SyntaxError: invalid syntax

@stefdoerr
Copy link
Contributor

I think the issue is simply that you put chain = " " meaning that the chain should be an empty space instead of an empty string.
Try with:

atom1_alpha_Dihang_residue_6 = {'name': "O3'", 'resid': 5, 'segid': 'RNAA', 'chain': ''}

Also make sure the segid is RNAA by doing print(mol.segid)

@vas2201
Copy link
Author

vas2201 commented Mar 22, 2023

The above issue is resolved in the Charmm force field.
However, when I switched to amber forcefield for RNA. I am encountering the following issues.
Looks like the default amber force field not recognizing the RNA. Is there any way that I can specify the amber force field for RNA?.

I define the default amber force fields as per documentation as follows.
topos_amber = amber.defaultTopo()
frcmods_amber = amber.defaultParam()
bmol_amber = amber.build(rna, topo=topos_amber, param=frcmods_amber, outdir='./build_amber')

Attached files below for your reference
https://drive.google.com/file/d/1HGo6KUalPOYgSrtD9vsBXov5FkvUNHIn/view?usp=share_link


BuildError Traceback (most recent call last)
Cell In[13], line 3
1 topos_amber = amber.defaultTopo()
2 frcmods_amber = amber.defaultParam()
----> 3 bmol_amber = amber.build(rna, topo=topos_amber, param=frcmods_amber, outdir='./build_amber')

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/builder/amber.py:469, in build(mol, ff, topo, param, prefix, outdir, caps, ionize, saltconc, saltanion, saltcation, disulfide, teleap, teleapimports, execute, atomtypes, offlibraries, gbsa, igb)
466 _detect_cofactors_ncaa_ptm(mol, param, topo)
468 if ionize:
--> 469 molbuilt = _build(
470 mol,
471 ff=ff,
472 topo=topo,
473 param=param,
474 prefix=prefix,
475 outdir=outdir,
476 disulfide=disulfide,
477 teleap=teleap,
478 teleapimports=teleapimports,
479 execute=execute,
480 atomtypes=atomtypes,
481 offlibraries=offlibraries,
482 gbsa=gbsa,
483 igb=igb,
484 )
485 shutil.move(
486 os.path.join(outdir, "structure.crd"),
487 os.path.join(outdir, "structure.noions.crd"),
488 )
489 shutil.move(
490 os.path.join(outdir, "structure.prmtop"),
491 os.path.join(outdir, "structure.noions.prmtop"),
492 )

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/builder/amber.py:728, in _build(mol, ff, topo, param, prefix, outdir, disulfide, teleap, teleapimports, execute, atomtypes, offlibraries, gbsa, igb)
726 os.chdir(currdir)
727 if errors:
--> 728 raise BuildError(
729 errors
730 + [f"Check {logpath} for further information on errors in building."]
731 )
732 logger.info("Finished building.")
734 if (
735 os.path.exists(os.path.join(outdir, "structure.crd"))
736 and os.path.getsize(os.path.join(outdir, "structure.crd")) != 0
737 and os.path.getsize(os.path.join(outdir, "structure.prmtop")) != 0
738 ):

BuildError: UnknownResidueError("Unknown residue(s) ['URA'] found in the input structure. You are either missing a topology definition for the residue or you need to rename it to the correct residue name")
MissingAtomTypeError('Missing atom type for [".R<DG5 1>.A<O2' 32>", ".R<DG5 1>.A<HO2' 33>", ".R<DG 2>.A<O2' 34>", ".R<DG 2>.A<HO2' 35>", ".R<DA 3>.A<O2' 33>", ".R<DA 3>.A<HO2' 34>", '.R<URA 4>.A<P 1>', '.R<URA 4>.A<O1P 2>', '.R<URA 4>.A<O2P 3>', ".R<URA 4>.A<O5' 4>", ".R<URA 4>.A<C5' 5>", ".R<URA 4>.A<C4' 6>", ".R<URA 4>.A<O4' 7>", ".R<URA 4>.A<C3' 8>", ".R<URA 4>.A<O3' 9>", ".R<URA 4>.A<C2' 10>", ".R<URA 4>.A<O2' 11>", ".R<URA 4>.A<C1' 12>", '.R<URA 4>.A<N1 13>', '.R<URA 4>.A<C2 14>', '.R<URA 4>.A<O2 15>', '.R<URA 4>.A<N3 16>', '.R<URA 4>.A<C4 17>', '.R<URA 4>.A<O4 18>', '.R<URA 4>.A<C5 19>', '.R<URA 4>.A<C6 20>', ".R<URA 4>.A<H5' 21>", ".R<URA 4>.A<H5'' 22>", ".R<URA 4>.A<H4' 23>", ".R<URA 4>.A<H3' 24>", ".R<URA 4>.A<H2' 25>", ".R<URA 4>.A<HO2' 26>", ".R<URA 4>.A<H1' 27>", '.R<URA 4>.A<H3 28>', '.R<URA 4>.A<H5 29>', '.R<URA 4>.A<H6 30>', ".R<DC 5>.A<O2' 31>", ".R<DC 5>.A<HO2' 32>", ".R<DA 6>.A<O2' 33>", ".R<DA 6>.A<HO2' 34>", ".R<DA 7>.A<O2' 33>", ".R<DA 7>.A<HO2' 34>", '.R<URA 8>.A<P 1>', '.R<URA 8>.A<O1P 2>', '.R<URA 8>.A<O2P 3>', ".R<URA 8>.A<O5' 4>", ".R<URA 8>.A<C5' 5>", ".R<URA 8>.A<C4' 6>", ".R<URA 8>.A<O4' 7>", ".R<URA 8>.A<C3' 8>", ".R<URA 8>.A<O3' 9>", ".R<URA 8>.A<C2' 10>", ".R<URA 8>.A<O2' 11>", ".R<URA 8>.A<C1' 12>", '.R<URA 8>.A<N1 13>', '.R<URA 8>.A<C2 14>', '.R<URA 8>.A<O2 15>', '.R<URA 8>.A<N3 16>', '.R<URA 8>.A<C4 17>', '.R<URA 8>.A<O4 18>', '.R<URA 8>.A<C5 19>', '.R<URA 8>.A<C6 20>', ".R<URA 8>.A<H5' 21>", ".R<URA 8>.A<H5'' 22>", ".R<URA 8>.A<H4' 23>", ".R<URA 8>.A<H3' 24>", ".R<URA 8>.A<H2' 25>", ".R<URA 8>.A<HO2' 26>", ".R<URA 8>.A<H1' 27>", '.R<URA 8>.A<H3 28>', '.R<URA 8>.A<H5 29>', '.R<URA 8>.A<H6 30>', ".R<DA 9>.A<O2' 33>", ".R<DA 9>.A<HO2' 34>", ".R<DG 10>.A<O2' 34>", ".R<DG 10>.A<HO2' 35>", ".R<DC 11>.A<O2' 31>", ".R<DC 11>.A<HO2' 32>", ".R<DA 12>.A<O2' 33>", ".R<DA 12>.A<HO2' 34>", ".R<DG 13>.A<O2' 34>", ".R<DG 13>.A<HO2' 35>", ".R<DG 14>.A<O2' 34>", ".R<DG 14>.A<HO2' 35>", '.R<URA 15>.A<P 1>', '.R<URA 15>.A<O1P 2>', '.R<URA 15>.A<O2P 3>', ".R<URA 15>.A<O5' 4>", ".R<URA 15>.A<C5' 5>", ".R<URA 15>.A<C4' 6>", ".R<URA 15>.A<O4' 7>", ".R<URA 15>.A<C3' 8>", ".R<URA 15>.A<O3' 9>", ".R<URA 15>.A<C2' 10>", ".R<URA 15>.A<O2' 11>", ".R<URA 15>.A<C1' 12>", '.R<URA 15>.A<N1 13>', '.R<URA 15>.A<C2 14>', '.R<URA 15>.A<O2 15>', '.R<URA 15>.A<N3 16>', '.R<URA 15>.A<C4 17>', '.R<URA 15>.A<O4 18>', '.R<URA 15>.A<C5 19>', '.R<URA 15>.A<C6 20>', ".R<URA 15>.A<H5' 21>", ".R<URA 15>.A<H5'' 22>", ".R<URA 15>.A<H4' 23>", ".R<URA 15>.A<H3' 24>", ".R<URA 15>.A<H2' 25>", ".R<URA 15>.A<HO2' 26>", ".R<URA 15>.A<H1' 27>", '.R<URA 15>.A<H3 28>', '.R<URA 15>.A<H5 29>', '.R<URA 15>.A<H6 30>', ".R<DG 16>.A<O2' 34>", ".R<DG 16>.A<HO2' 35>", '.R<URA 17>.A<P 1>', '.R<URA 17>.A<O1P 2>', '.R<URA 17>.A<O2P 3>', ".R<URA 17>.A<O5' 4>", ".R<URA 17>.A<C5' 5>", ".R<URA 17>.A<C4' 6>", ".R<URA 17>.A<O4' 7>", ".R<URA 17>.A<C3' 8>", ".R<URA 17>.A<O3' 9>", ".R<URA 17>.A<C2' 10>", ".R<URA 17>.A<O2' 11>", ".R<URA 17>.A<C1' 12>", '.R<URA 17>.A<N1 13>', '.R<URA 17>.A<C2 14>', '.R<URA 17>.A<O2 15>', '.R<URA 17>.A<N3 16>', '.R<URA 17>.A<C4 17>', '.R<URA 17>.A<O4 18>', '.R<URA 17>.A<C5 19>', '.R<URA 17>.A<C6 20>', ".R<URA 17>.A<H5' 21>", ".R<URA 17>.A<H5'' 22>", ".R<URA 17>.A<H4' 23>", ".R<URA 17>.A<H3' 24>", ".R<URA 17>.A<H2' 25>", ".R<URA 17>.A<HO2' 26>", ".R<URA 17>.A<H1' 27>", '.R<URA 17>.A<H3 28>', '.R<URA 17>.A<H5 29>', '.R<URA 17>.A<H6 30>', ".R<DG 18>.A<O2' 34>", ".R<DG 18>.A<HO2' 35>", ".R<DG 19>.A<O2' 34>", ".R<DG 19>.A<HO2' 35>", ".R<DC 20>.A<O2' 31>", ".R<DC 20>.A<HO2' 32>", ".R<DA 21>.A<O2' 33>", ".R<DA 21>.A<HO2' 34>", ".R<DC 22>.A<O2' 31>", ".R<DC 22>.A<HO2' 32>", ".R<DA 23>.A<O2' 33>", ".R<DA 23>.A<HO2' 34>", ".R<DC 24>.A<O2' 31>", ".R<DC 24>.A<HO2' 32>", ".R<DC 25>.A<O2' 31>", ".R<DC 25>.A<HO2' 32>", ".R<DA 26>.A<O2' 33>", ".R<DA 26>.A<HO2' 34>", ".R<DG 27>.A<O2' 34>", ".R<DG 27>.A<HO2' 35>", '.R<URA 28>.A<P 1>', '.R<URA 28>.A<O1P 2>', '.R<URA 28>.A<O2P 3>', ".R<URA 28>.A<O5' 4>", ".R<URA 28>.A<C5' 5>", ".R<URA 28>.A<C4' 6>", ".R<URA 28>.A<O4' 7>", ".R<URA 28>.A<C3' 8>", ".R<URA 28>.A<O3' 9>", ".R<URA 28>.A<C2' 10>", ".R<URA 28>.A<O2' 11>", ".R<URA 28>.A<C1' 12>", '.R<URA 28>.A<N1 13>', '.R<URA 28>.A<C2 14>', '.R<URA 28>.A<O2 15>', '.R<URA 28>.A<N3 16>', '.R<URA 28>.A<C4 17>', '.R<URA 28>.A<O4 18>', '.R<URA 28>.A<C5 19>', '.R<URA 28>.A<C6 20>', ".R<URA 28>.A<H5' 21>", ".R<URA 28>.A<H5'' 22>", ".R<URA 28>.A<H4' 23>", ".R<URA 28>.A<H3' 24>", ".R<URA 28>.A<H2' 25>", ".R<URA 28>.A<HO2' 26>", ".R<URA 28>.A<H1' 27>", '.R<URA 28>.A<H3 28>', '.R<URA 28>.A<H5 29>', '.R<URA 28>.A<H6 30>', ".R<DC 29>.A<O2' 31>", ".R<DC 29>.A<HO2' 32>", ".R<DA 30>.A<O2' 33>", ".R<DA 30>.A<HO2' 34>", '.R<URA 31>.A<P 1>', '.R<URA 31>.A<O1P 2>', '.R<URA 31>.A<O2P 3>', ".R<URA 31>.A<O5' 4>", ".R<URA 31>.A<C5' 5>", ".R<URA 31>.A<C4' 6>", ".R<URA 31>.A<O4' 7>", ".R<URA 31>.A<C3' 8>", ".R<URA 31>.A<O3' 9>", ".R<URA 31>.A<C2' 10>", ".R<URA 31>.A<O2' 11>", ".R<URA 31>.A<C1' 12>", '.R<URA 31>.A<N1 13>', '.R<URA 31>.A<C2 14>', '.R<URA 31>.A<O2 15>', '.R<URA 31>.A<N3 16>', '.R<URA 31>.A<C4 17>', '.R<URA 31>.A<O4 18>', '.R<URA 31>.A<C5 19>', '.R<URA 31>.A<C6 20>', ".R<URA 31>.A<H5' 21>", ".R<URA 31>.A<H5'' 22>", ".R<URA 31>.A<H4' 23>", ".R<URA 31>.A<H3' 24>", ".R<URA 31>.A<H2' 25>", ".R<URA 31>.A<HO2' 26>", ".R<URA 31>.A<H1' 27>", '.R<URA 31>.A<H3 28>', '.R<URA 31>.A<H5 29>', '.R<URA 31>.A<H6 30>', ".R<DA 32>.A<O2' 33>", ".R<DA 32>.A<HO2' 34>", ".R<DC 33>.A<O2' 31>", ".R<DC 33>.A<HO2' 32>", ".R<DC 34>.A<O2' 31>", ".R<DC 34>.A<HO2' 32>", '.R<URA 35>.A<P 1>', '.R<URA 35>.A<O1P 2>', '.R<URA 35>.A<O2P 3>', ".R<URA 35>.A<O5' 4>", ".R<URA 35>.A<C5' 5>", ".R<URA 35>.A<C4' 6>", ".R<URA 35>.A<O4' 7>", ".R<URA 35>.A<C3' 8>", ".R<URA 35>.A<O3' 9>", ".R<URA 35>.A<C2' 10>", ".R<URA 35>.A<O2' 11>", ".R<URA 35>.A<C1' 12>", '.R<URA 35>.A<N1 13>', '.R<URA 35>.A<C2 14>', '.R<URA 35>.A<O2 15>', '.R<URA 35>.A<N3 16>', '.R<URA 35>.A<C4 17>', '.R<URA 35>.A<O4 18>', '.R<URA 35>.A<C5 19>', '.R<URA 35>.A<C6 20>', ".R<URA 35>.A<H5' 21>", ".R<URA 35>.A<H5'' 22>", ".R<URA 35>.A<H4' 23>", ".R<URA 35>.A<H3' 24>", ".R<URA 35>.A<H2' 25>", ".R<URA 35>.A<HO2' 26>", ".R<URA 35>.A<H1' 27>", '.R<URA 35>.A<H3 28>', '.R<URA 35>.A<H5 29>', '.R<URA 35>.A<H6 30>', '.R<URA 36>.A<P 1>', '.R<URA 36>.A<O1P 2>', '.R<URA 36>.A<O2P 3>', ".R<URA 36>.A<O5' 4>", ".R<URA 36>.A<C5' 5>", ".R<URA 36>.A<C4' 6>", ".R<URA 36>.A<O4' 7>", ".R<URA 36>.A<C3' 8>", ".R<URA 36>.A<O3' 9>", ".R<URA 36>.A<C2' 10>", ".R<URA 36>.A<O2' 11>", ".R<URA 36>.A<C1' 12>", '.R<URA 36>.A<N1 13>', '.R<URA 36>.A<C2 14>', '.R<URA 36>.A<O2 15>', '.R<URA 36>.A<N3 16>', '.R<URA 36>.A<C4 17>', '.R<URA 36>.A<O4 18>', '.R<URA 36>.A<C5 19>', '.R<URA 36>.A<C6 20>', ".R<URA 36>.A<H5' 21>", ".R<URA 36>.A<H5'' 22>", ".R<URA 36>.A<H4' 23>", ".R<URA 36>.A<H3' 24>", ".R<URA 36>.A<H2' 25>", ".R<URA 36>.A<HO2' 26>", ".R<URA 36>.A<H1' 27>", '.R<URA 36>.A<H3 28>', '.R<URA 36>.A<H5 29>', '.R<URA 36>.A<H6 30>', ".R<DG 37>.A<O2' 34>", ".R<DG 37>.A<HO2' 35>", ".R<DA 38>.A<O2' 33>", ".R<DA 38>.A<HO2' 34>", '.R<URA 39>.A<P 1>', '.R<URA 39>.A<O1P 2>', '.R<URA 39>.A<O2P 3>', ".R<URA 39>.A<O5' 4>", ".R<URA 39>.A<C5' 5>", ".R<URA 39>.A<C4' 6>", ".R<URA 39>.A<O4' 7>", ".R<URA 39>.A<C3' 8>", ".R<URA 39>.A<O3' 9>", ".R<URA 39>.A<C2' 10>", ".R<URA 39>.A<O2' 11>", ".R<URA 39>.A<C1' 12>", '.R<URA 39>.A<N1 13>', '.R<URA 39>.A<C2 14>', '.R<URA 39>.A<O2 15>', '.R<URA 39>.A<N3 16>', '.R<URA 39>.A<C4 17>', '.R<URA 39>.A<O4 18>', '.R<URA 39>.A<C5 19>', '.R<URA 39>.A<C6 20>', ".R<URA 39>.A<H5' 21>", ".R<URA 39>.A<H5'' 22>", ".R<URA 39>.A<H4' 23>", ".R<URA 39>.A<H3' 24>", ".R<URA 39>.A<H2' 25>", ".R<URA 39>.A<HO2' 26>", ".R<URA 39>.A<H1' 27>", '.R<URA 39>.A<H3 28>', '.R<URA 39>.A<H5 29>', '.R<URA 39>.A<H6 30>', ".R<DC 40>.A<O2' 31>", ".R<DC 40>.A<HO2' 32>", ".R<DC3 41>.A<O2' 32>", ".R<DC3 41>.A<HO2' 33>"]')
Check /home/nmrbox/spenumutchu/Desktop/Projects/aufr12_project/SL2_drug/build_amber/log.txt for further information on errors in building.

from htmd.ui import *

@stefdoerr
Copy link
Contributor

Can you send me the initial file you used? Normally you should be able to build RNA fine in AMBER.

@vas2201
Copy link
Author

vas2201 commented Mar 23, 2023

I uploaded the file that i used with charmm force field for RNA. Running fine without issues.

https://drive.google.com/file/d/1R3DajOODBxT621giw9RVkaNBgyMOtSWV/view?usp=share_link

@stefdoerr
Copy link
Contributor

No sorry I meant the starting PDB file. Because I want to test the pipeline from the start

@vas2201
Copy link
Author

vas2201 commented Mar 23, 2023

@stefdoerr
Copy link
Contributor

There is something non-standard about your atom or residue names that pdb2pqr does not like and it doesn't recognize the nucleic acids as such. If I use your file I get:

ValueError: No biomolecule heavy atoms found and no ligand present. Unable to proceed.  You may also see this message if PDB2PQR does not have parameters for any residue in your biomolecule.

If I used the 5v16 directly from the RCSB database it works correctly

from htmd.builder.amber import build
from htmd.builder.solvate import solvate
from moleculekit.tools.preparation import systemPrepare
from moleculekit.molecule import Molecule

mol = Molecule("5v16")
pmol = systemPrepare(mol)
smol = solvate(pmol, pad=10)
bmol = build(smol, ionize=True, outdir="/tmp/test")

I don't recommend solvating like that though (you should solvate as a cubic box) but just wanted to demonstrate that it works fine if you start from the RCSB file

@vas2201
Copy link
Author

vas2201 commented Mar 31, 2023

I tried using adaptive sampling as per your suggestions with the Amber force field and user-defined dihedrals for RNA. While I did not encounter any issues with the Charmm force field and adaptive sampling was successful, I am having issues with adaptive sampling using the Amber force field for RNA (after first epoch). I tried tweaking the code for segid and chain in multiple ways to resolve the issues, but encountered issues such as atoms not being selected or detected. If you have time, it would be great if you could review the code. The Jupyter notebook can be found at this link: https://drive.google.com/file/d/1WEOASnxw62U44T-aSfD1SH70aoR1ISCE/view?usp=share_link

Error : 
2023-03-30 01:28:08,531 - jobqueues.localqueue - INFO - Running /home/nmrbox/spenumutchu/Desktop/Projects/aufr12_project/SL2_WT/htmd_amber/input/e1s9_md_50ns_1 on device 1
2023-03-30 01:28:30,000 - htmd.adaptive.adaptive - INFO - Processing epoch 1
2023-03-30 01:28:30,000 - htmd.adaptive.adaptive - INFO - Retrieving simulations.
2023-03-30 01:28:30,001 - htmd.adaptive.adaptive - INFO - 2 simulations in progress
2023-03-30 01:28:30,001 - htmd.adaptive.adaptiverun - INFO - Postprocessing new data
Creating simlist: 100%|███████████████████████████| 9/9 [00:00<00:00, 63.93it/s]
2023-03-30 01:28:32,094 - htmd.simlist - WARNING - Filtering was not able to write ./filtered/filtered.prmtop due to error: Molecule cannot write files with "prmtop" extension yet. If you need such support please notify us on the github moleculekit issue tracker.
Filtering trajectories: 100%|█████████████████████| 8/8 [00:09<00:00,  1.25s/it]

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[20], line 2
      1 get_ipython().run_line_magic('cd', '/home/nmrbox/spenumutchu/Desktop/Projects/aufr12_project/SL2_WT/htmd_amber')
----> 2 ad.run()

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/adaptive/adaptive.py:196, in AdaptiveBase.run(self)
    194 # If currently running simulations are lower than nmin start new ones to reach nmax number of sims
    195 if self._running <= self.nmin and epoch < self.nepochs:
--> 196     flag = self._algorithm()
    197     if flag is False:
    198         self._unsetLock()

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/adaptive/adaptiverun.py:208, in AdaptiveMD._algorithm(self)
    207 def _algorithm(self):
--> 208     data = self._getData(self._getSimlist())
    209     if not self._checkNFrames(data):
    210         return False

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/adaptive/adaptiverun.py:246, in AdaptiveMD._getData(self, sims)
    241 # if self.contactsym is not None:
    242 #    contactSymmetry(data, self.contactsym)
    244 if self.ticadim > 0:
    245     # tica = TICA(metr, int(max(2, np.ceil(self.ticalag))))  # gianni: without project it was tooooo slow
--> 246     data = metr.project()
    247     data.dropTraj()  # Drop before TICA to avoid broken trajectories
    248     # 1 < ticalag < (trajLen / 2)

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/htmd/projections/metric.py:159, in Metric.project(self, njobs)
    157     for proj in self.projectionlist:
    158         if isinstance(proj, Projection):
--> 159             proj._setCache(uqMol)
    160 else:
    161     logger.warning(
    162         "Cannot calculate description of dimensions due to different topology files for each trajectory."
    163     )

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/moleculekit/projections/projection.py:31, in Projection._setCache(self, mol)
     30 def _setCache(self, mol):
---> 31     resdict = self._calculateMolProp(mol)
     32     self._cache.update(resdict)

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/moleculekit/projections/metricdihedral.py:839, in MetricDihedral._calculateMolProp(self, mol, props)
    836 else:
    837     dihedrals = ensurelist(self._dihedrals)
--> 839 res["dihedrals"] = Dihedral.dihedralsToIndexes(mol, dihedrals, protsel)
    840 return res

File ~/anaconda3/envs/100htmd/lib/python3.9/site-packages/moleculekit/projections/metricdihedral.py:136, in Dihedral.dihedralsToIndexes(mol, dihedrals, sel)
    134     atomsel = atomsel & selatoms
    135     if np.sum(atomsel) != 1:
--> 136         raise RuntimeError(
    137             "Expected one atom from atomselection {}. Got {} instead.".format(
    138                 a, np.sum(atomsel)
    139             )
    140         )
    141     idx.append(np.where(atomsel)[0][0])
    142 indexes.append(idx)

RuntimeError: Expected one atom from atomselection {'name': "O3'", 'resid': 5, 'chain': '', 'insertion': '', 'segid': ''}. Got 0 instead.

@stefdoerr
Copy link
Contributor

Can you share with me these folders please so I can take a look?

ad.generatorspath = './generators'
ad.inputpath = './input'
ad.datapath = './data'
ad.filteredpath = './filtered'

@vas2201
Copy link
Author

vas2201 commented Apr 4, 2023

Thanks for your reply. I attached zip folder as drive link and contain requested folders.
https://drive.google.com/file/d/10qc_C5LH-ZkUXjA9WHbXkLK0knQv7S_G/view?usp=share_link
Regards
Vas

@stefdoerr
Copy link
Contributor

After filtering the Molecule has a segid "0"

In [16]: mol.segid
Out[16]: array(['0', '0', '0', ..., '0', '0', '0'], dtype=object)

So the error it gives

RuntimeError: Expected one atom from atomselection {'name': "O3'", 'resid': 5, 'chain': '', 'insertion': '', 'segid': ''}. Got 0 instead.

is correct because segid is not "".
Fix your dihedral selections to have "segid": "0"

@vas2201
Copy link
Author

vas2201 commented Apr 5, 2023

I have identified the problem. Upon building the build_amber, the segid was no longer present . That is causing conflict while defining the dihedrals in metric.

To address this, I utilized the following command to define the segid , which successfully resolved the issue.

bmol.set('segid', '0', sel='nucleic')

Thanks for your help on this matter.
Regards
Vas

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants