search

Acellera / htmd

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

calculating interaction energy in ACEMD #1006

Closed zainbsa closed 2 years ago

zainbsa commented 3 years ago

I am new at working with ACEMD and still trying to figure things out, I will like to know if it is possible to calculate the interaction energies of ligands with protein and also plot/analyze the contacts of the ligand with a protein?

stefdoerr commented 3 years ago

Hi, with ACEMD you cannot get specific interaction energies out. It is however possible with ffevaluate https://software.acellera.com/docs/latest/htmd/tutorials/FFEvaluate.html

To calculate and plot protein ligand contacts you can do it with moleculekit MetricDistance class: https://software.acellera.com/docs/latest/moleculekit/moleculekit.projections.metricdistance.html You pass it a Molecule object with your trajectory and two selections for which to calculate distances and it returns those distances

zainbsa commented 3 years ago

Ok! Thank you.I will try these out.

zainbsa commented 3 years ago

Hi, I followed the sample at https://software.acellera.com/docs/latest/htmd/tutorials/FFEvaluate.html to evaluate the interaction energies and setup my job script, but I have the topology extension as .top, and while running this I got the error below; RuntimeError: Extension of file /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/structure.top not recognized. Report issue on HTMD issue tracker.

Does it mean it only recognize the .prmtop extension as in the sample? Please note I had to convert my topology from charmm to amber using the below; import parmed as pmd pmd.load_file('structure.psf').save('structure.top')

I also tried to change the extension here to .prmtop but it keeps refusing with this error; AttributeError: 'NoneType' object has no attribute 'idx'

stefdoerr commented 3 years ago

Hi, unfortunately you will have to find a way to convert to prmtop. Maybe antechamber can do it, otherwise you might have to rebuild the system in AMBER. Sorry I cannot help you further but I'm currently out on holidays.

On Fri, Aug 13, 2021, 16:05 zainbsa @.***> wrote:

Hi, I followed the sample at https://software.acellera.com/docs/latest/htmd/tutorials/FFEvaluate.html to evaluate the interaction energies and setup my job script, but I have the topology extension as .top, and while running this I got the error below; RuntimeError: Extension of file /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/structure.top not recognized. Report issue on HTMD issue tracker.

Does it mean it only recognize the .prmtop extension as in the sample? Please note I had to convert my topology from charmm to amber using the below; import parmed as pmd pmd.load_file('structure.psf').save('structure.top')

I also tried to change the extension here to .prmtop but it keeps refusing with this error; AttributeError: 'NoneType' object has no attribute 'idx'

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/Acellera/htmd/issues/1006#issuecomment-898443740, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RLASZ2S5L4NVLXQZ2AZ3T4UKARANCNFSM5B6GTNNQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

zainbsa commented 3 years ago

Hi, I suspect the conversion I made from .psf file to .top using; pmd.load_file('structure.psf').save('structure.top') as initially stated is to a gromacs topology file and that is why it doesn't recognize the extension file while trying out the sample here; https://software.acellera.com/docs/latest/htmd/tutorials/FFEvaluate.html

However, I tried to covert the topology file in gromacs to amber using; gmx_top = pmd.load_file('structure.top', xyz= 'filtered.gro') gmx_top.save('name.top', format='amber') gmx_top.save('name.crd', format='rst7')

and then got this error; MoleculeError: Cannot bond two virtual sites/extra points together. couldnt figure the conversion with antechamber

zainbsa commented 3 years ago

Hello,

As I understand now using the FFEvaluate the loadParameters function is suppose to read both AMBER and CHARMM forcefields but while using the .psf file generated from CHARMM I keep getting an error that the extension is not recognized. RuntimeError: Extension of file /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/structure.psf not recognized. Report issue on HTMD issue tracker.

stefdoerr commented 3 years ago

This is weird. Can you post the code which throws the error and the full traceback which it prints? Thanks!

zainbsa commented 3 years ago

Here it is; 

from ffevaluation.ffevaluate import FFEvaluate, loadParameters, viewForces
from moleculekit.molecule import Molecule
from htmd.home import home
from os.path import join

datadir = home(dataDir='abeta_sergio')

# Load the parameters
prm = loadParameters(join(datadir, 'structure.psf'))
# Load the topology
mol = Molecule(join(datadir, 'structure.psf'))

RuntimeError                              Traceback (most recent call last)
<ipython-input-18-2d29d01da3ab> in <module>
      9 
     10 # Load the parameters
---> 11 prm = loadParameters(join(datadir, 'structure.psf'))
     12 # Load the topology
     13 mol = Molecule(join(datadir, 'structure.psf'))

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in loadParameters(fname)
     70 
     71     if prm is None:
---> 72         raise RuntimeError('Extension of file {} not recognized. Report issue on HTMD issue tracker.'.format(fname))
     73 
     74     return prm

RuntimeError: Extension of file /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/structure.psf not recognized. Report issue on HTMD issue tracker.
stefdoerr commented 3 years ago

PSF files don't contain parameters, they only contain topology. You also need a prm which the loadParameters function expects.

zainbsa commented 3 years ago

I see.However, using a parameter file; I still got the same error

# Load the parameters
prm = loadParameters(join(datadir, 'par_all36m_prot.prm'))
# Load the topology
mol = Molecule(join(datadir, 'structure.psf'))

RuntimeError                              Traceback (most recent call last)
<ipython-input-21-5184fbc10d0d> in <module>
      8 
      9 # Load the parameters
---> 10 prm = loadParameters(join(datadir, 'par_all36m_prot.prm'))
     11 # Load the topology
     12 mol = Molecule(join(datadir, 'structure.psf'))

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in loadParameters(fname)
     70 
     71     if prm is None:
---> 72         raise RuntimeError('Extension of file {} not recognized. Report issue on HTMD issue tracker.'.format(fname))
     73 
     74     return prm

RuntimeError: Extension of file /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/par_all36m_prot.prm not recognized. Report issue on HTMD issue tracker.
stefdoerr commented 3 years ago

Did it print anything else before that complaining that it could not read the file correctly?

stefdoerr commented 3 years ago

Because the function definitely can read prm files https://github.com/Acellera/ffevaluation/blob/master/ffevaluation/ffevaluate.py#L33-L74

zainbsa commented 3 years ago

yes... Failed to read /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/tmp_gauss-min.prm as CHARMM parameters. Attempting with AMBER prmtop reader Failed to read /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/tmp_gauss-min.prm due to errors [Errno 2] No such file or directory: '/software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/tmp_gauss-min.prm' Could not open /software/software/HTMD/1.23.5/lib/python3.6/site-packages/htmd/data/abeta_sergio/tmp_gauss-min.prm for reading

stefdoerr commented 3 years ago

Then it's probably not a valid prm file. Try this and see what it gives

import parmed
prm = parmed.charmm.CharmmParameterSet("myfile.prm")
zainbsa commented 3 years ago

Using the above command line it actually reads the file, I also had to load my topology this way; mol = parmed.charmm.CharmmParameterSet ('0.top_all36_prot.rtf') (P.S it doesn't read the .psf file here)

However, for the ffevaluation (I got the below error, could you please help me rewrite this script?)

# Create a FFEvaluate object
ffev = FFEvaluate(mol, prm)

AttributeError                            Traceback (most recent call last)
<ipython-input-36-f61058ac0bbd> in <module>
     15 
     16 # Create a FFEvaluate object
---> 17 ffev = FFEvaluate(mol, prm)
     18 
     19 # Load some coordinates

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in __init__(self, mol, prm, betweensets, cutoff, rfa, solventDielectric)
    138         >>> energies, forces, atmnrg = ffev.calculate(mol.coords)
    139         """
--> 140         mol = mol.copy()
    141         setA, setB = calculateSets(mol, betweensets)
    142 

AttributeError: 'CharmmParameterSet' object has no attribute 'copy'
stefdoerr commented 3 years ago 
mol = Molecule("structure.psf")
prm = parmed.charmm.CharmmParameterSet("myfile.prm")
ffev = FFEvaluate(mol, prm)
zainbsa commented 3 years ago

I am sorry but it doesn't get better

KeyError                                  Traceback (most recent call last)
<ipython-input-37-3f0e6ce59fd6> in <module>
     15 
     16 # Create a FFEvaluate object
---> 17 ffev = FFEvaluate(mol, prm)
     18 
     19 # Load some coordinates

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in __init__(self, mol, prm, betweensets, cutoff, rfa, solventDielectric)
    141         setA, setB = calculateSets(mol, betweensets)
    142 
--> 143         args = list(init(mol, prm))
    144         args.append(setA)
    145         args.append(setB)

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in init(mol, prm)
    285     epsilon14 = np.zeros(len(uqtypes), dtype=np.float32)
    286     for i, type in enumerate(uqtypes):
--> 287         sigma[i] = prm.atom_types[type].sigma
    288         epsilon[i] = prm.atom_types[type].epsilon
    289         sigma14[i] = prm.atom_types[type].sigma_14

KeyError: 'CLA'
stefdoerr commented 3 years ago

There are missing some parameters in your prm file. Specifically those for residue names CLA which is the chloride ion. Try running

from ffevaluation.test_ffevaluate import fixParameters
fixParameters("myprmfile.prm", "myprmfile_fixed.prm")

mol = Molecule("structure.psf")
prm = parmed.charmm.CharmmParameterSet("myprmfile_fixed.prm")
ffev = FFEvaluate(mol, prm)

If this doesn't work then you might need to copy those parameters into the prm file

zainbsa commented 3 years ago

Hello, I loaded all the parameter files and I do not have the above error anymore but came up with another error. Is the complain still about the parameters? I am so sorry I just really need to make this work.

prm = parmed.charmm.CharmmParameterSet('par_all36m_prot.prm', 'tmp_gauss-min_charmm.rtf', 'toppar_water_ions_test.str', 'tmp_gauss-min_charmm.par', 'toppar_water_ions.str', 'top_water_ions.rtf')

mol = Molecule('structure.psf')

# Create a FFEvaluate object
ffev = FFEvaluate(mol, prm)

# Load some coordinates
mol.read(join(datadir,'structure.pdb'))

RuntimeError                              Traceback (most recent call last)
<ipython-input-16-32a9d048c7d7> in <module>
     18 
     19 # Create a FFEvaluate object
---> 20 ffev = FFEvaluate(mol, prm)
     21 
     22 # Load some coordinates

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in __init__(self, mol, prm, betweensets, cutoff, rfa, solventDielectric)
    141         setA, setB = calculateSets(mol, betweensets)
    142 
--> 143         args = list(init(mol, prm))
    144         args.append(setA)
    145         args.append(setB)

/software/software/HTMD/1.23.5/lib/python3.6/site-packages/ffevaluation/ffevaluate.py in init(mol, prm)
    340             dihparam = prm.dihedral_types[ty[::-1]]
    341         else:
--> 342             raise RuntimeError('Could not find type {} in dihedral_types'.format(ty))
    343         i, j = sorted([dihed[0], dihed[3]])
    344         s14_atom_list[i].append(j)

RuntimeError: Could not find type ('NH1', 'CT1', 'CT3', 'HA3') in dihedral_types
stefdoerr commented 3 years ago

Is the gauss min file a ligand? If yes can you check their prm files if there are atom types defined called NH1 CT1 CT3 HA3? It seems to be missing the energy terms for the above dihedral angle.

zainbsa commented 3 years ago

Yea It is a ligand and I check the prm files it doesnt have the atom defined called NH1 CT1 CT3 ...How do I resolve this?

BONDS HNRP NRP 443.528 1.0280 HCMM CR 342.991 1.0930 NRP CR 276.638 1.4800 CR CR 306.432 1.5080 O2CM SO2 773.493 1.4500 CR SO2 234.466 1.7720

ANGLES CR NRP HNRP 41.452 111.2060 HNRP NRP HNRP 41.596 107.7870 NRP CR CR 84.848 106.4930 NRP CR HCMM 62.754 106.2240 CR CR HCMM 45.770 110.5490 HCMM CR HCMM 37.134 108.8360 CR CR CR 61.243 109.6080 O2CM SO2 O2CM 112.914 120.9240 O2CM SO2 CR 104.062 107.0660 CR CR SO2 78.658 109.3150 SO2 CR HCMM 47.713 106.8550 DIHEDRALS NRP CR CR CR -0.324 1 0.00 NRP CR CR CR 0.275 2 180.00 NRP CR CR CR 0.295 3 0.00 NRP CR CR HCMM 0.346 1 0.00 NRP CR CR HCMM -0.265 2 180.00 NRP CR CR HCMM 0.139 3 0.00 CR CR CR SO2 0.150 3 0.00 CR CR CR HCMM 0.320 1 0.00 CR CR CR HCMM -0.315 2 180.00 CR CR CR HCMM 0.132 3 0.00 CR CR NRP HNRP 0.093 3 0.00 CR CR SO2 O2CM 0.050 3 0.00 SO2 CR CR HCMM 0.150 3 0.00 O2CM SO2 CR HCMM 0.292 2 180.00 O2CM SO2 CR HCMM 0.194 3 0.00 HNRP NRP CR HCMM 0.130 3 0.00 HCMM CR CR HCMM 0.142 1 0.00 HCMM CR CR HCMM -0.693 2 180.00 HCMM CR CR HCMM 0.157 3 0.00

IMPROPER NRP HNRP CR HNRP 0.000 0 0.00 CR CR NRP HCMM 0.000 0 0.00 CR HCMM NRP HCMM 0.000 0 0.00 CR CR CR HCMM 0.000 0 0.00 CR HCMM CR HCMM 0.000 0 0.00 CR SO2 CR HCMM 0.000 0 0.00 SO2 O2CM CR O2CM 0.000 0 0.00

NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch - cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5 NRP 0.000000 -0.200000 1.850000
CR 0.000000 -0.055000 2.175000 0.000000 -0.010000 1.900000
SO2 0.000000 -0.470000 2.100000
O2CM 0.000000 -0.120000 1.700000
HNRP 0.000000 -0.046000 0.224500
HCMM 0.000000 -0.022000 1.320000

stefdoerr commented 3 years ago

It's really hard to resolve such kind of issues. You would need to check which dihedral this corresponds to and then investigate why parmed cannot find parameters for it in your files. It's possible that it's an issue of parmed, or you might indeed be missing the dihedral parameters. If you want send me the files and I can take a look next week if I find some time.

zainbsa commented 3 years ago

I have a silly question before I send you the files when I tried to run this line; mol.read(join(datadir, 'structure.pdb')) I got a SyntaxError: invalid syntax ???

stefdoerr commented 3 years ago

The error is probably somewhere else. This command looks ok

On Thu, Aug 26, 2021, 12:58 zainbsa @.***> wrote:

I have a silly question before I send you the files when I tried to run this line; mol.read(join(datadir, 'structure.pdb')) I got a SyntaxError: invalid syntax ???

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/Acellera/htmd/issues/1006#issuecomment-906264370, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RLAS7QXAAEHA6UEPIQHDT6YF3TANCNFSM5B6GTNNQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

zainbsa commented 3 years ago

Kindly find the link for my files. The filtered. is after striping the waters and ions, while the parameters used for preparing the complex are tmp_gauss., top_water, par_all36m. Please let me know if I am missing anything. Thanks https://filesender.cesnet.cz/?s=download&token=86152803-7326-4ea2-9d22-6b2c19e8ffc2