htmd
htmd copied to clipboard
Acellera
RuntimeError (adaptive MD for RNA): when using a dihedral with an atom name containing a single quote
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.
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
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)
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 *
Can you send me the initial file you used? Normally you should be able to build RNA fine in AMBER.
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
No sorry I meant the starting PDB file. Because I want to test the pipeline from the start
attached pdb file here ! https://drive.google.com/file/d/1YWgrH4AyTHMP323EzXn5Ld2Gb0vxh5fP/view?usp=share_link
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
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.
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'
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
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"
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