H bond with dump LAMMPS trajectory

Jose Martín-Roca

unread,

Oct 24, 2022, 5:53:42 PM10/24/22

to MDnalysis discussion

Hi,

I did simulations of water (using tip4p2005 model) with LAMMPS long ago, but I used "dump custom" command to save my trajectory. I have been reading the documentation and mail list, and at the end I was able to load my trajectories using the following code:

import MDAnalysis as mda
u = mda.Universe("lammps.data",topology_format="LAMMPSDUMP", dt=0.005, atom_style="id type xs ys zs")

Then I have tried to compute the H bonds using:

from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis

hbonds = HydrogenBondAnalysis(u)

hbonds.run(verbose=True)

The problem is that I have not saved the charges of my atoms, and I get the following error

NoDataError: This Universe does not contain charge information

Since the molecules are just made of H and O, I would like to somehow set the value of the charge to this atoms to skip this error easily. How could I do this?

I have though that creating groups could be a possibility

Hidrogenos = u.select_atoms("type 1")
Oxigenos = u.select_atoms("type 2")

But I cannot find a way to add the charges to the Universe variable using this groups

Just in case someone need it, the trajectory has the following form:

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
648
ITEM: BOX BOUNDS pp pp pp
0.0000000000000000e+00 1.6649639999999998e+01
0.0000000000000000e+00 1.6649639999999998e+01
0.0000000000000000e+00 1.6649639999999998e+01
ITEM: ATOMS id type xs ys zs
153 1 0.085948 0.148918 0.0426951
353 1 0.290313 0.114785 0.128039
354 1 0.203023 0.139726 0.122935
216 1 0.197995 0.311325 0.139853
351 1 0.140253 0.0286082 0.13541
215 1 0.206278 0.399694 0.120106
349 2 0.0878958 0.0356515 0.158088
152 1 0.105248 0.234876 0.0201904
300 1 0.13702 0.351859 0.00795395

......

Thanks,

José

Oliver Beckstein

unread,

Oct 24, 2022, 7:48:45 PM10/24/22

to mdnalysis-discussion

Hi José,

You can add a charges attribute to the topology

u.add_TopologyAttr("charges”)

and then add the charges for the types

Hidrogenos.charges = 0.41

Oxigenos.charges = -0.82

Best,

Oliver

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/966ef826-b96c-4d1f-b4a4-2bcaebc988ean%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

Jose Martín-Roca

unread,

Oct 25, 2022, 2:34:49 PM10/25/22

to MDnalysis discussion

Hi Oliver,

Thank you for your answer, it works and it show an error with names and residues. I completed my code as you suggest and also added names and residues (see below). Maybe it is a trivial problem, but I have never work before with atomistic simulations. My problem now is that I don't know how to add the residues names to the universe variable, becase hbonds.run() still gives an error

SelectionError: Selection failed: 'Unexpected token 'and''

I guess this is due to the fact that residues is empty, isn't it?

ciao

José

-----------------------------------------------------------------------------------------------------------------------------------------------------------------

from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis

import MDAnalysis as mda
u = mda.Universe("lammps.data",topology_format="LAMMPSDUMP", dt=0.005, atom_style="id type xs ys zs")

Hidrogenos = u.select_atoms("type 1")
Oxigenos = u.select_atoms("type 2")

###############################
########### Masses #########
###############################
u.add_TopologyAttr("masses")
Hidrogenos.masses = 1.00794
Oxigenos.masses = 15.9994

###############################
########### Charges #########
###############################
u.add_TopologyAttr("charges")

Hidrogenos.charges = 0.41
Oxigenos.charges = -0.82

###############################
########### Names #########
###############################
nameH=[ ]
for ii in range(0,len(Hidrogenos)):
nameH.append('H')
nameO=[ ]
for ii in range(0,len(Oxigenos)):
nameO.append('O')

u.add_TopologyAttr("names")
Hidrogenos.names = nameH
Oxigenos.names = nameO

###############################
########### Residues ########
###############################
u.add_TopologyAttr("resnames")

###############################
####### Compute H bonds #######
###############################

hbonds = HydrogenBondAnalysis(u)

hbonds.run(verbose=True)

Oliver Beckstein

unread,

Oct 25, 2022, 6:45:18 PM10/25/22

to mdnalysis-discussion

HI José,

Not sure about where hbonds is failing.

You can easily add names like this:

water = u.select_atoms(“type 1 2”)

water.residues.resname = “HOH” # or whatever you want to name it

resnames should be assigned at the residues level. (You don’t need to make an array of names if you assign the same name to all.)

The hbonds user guide docs https://userguide.mdanalysis.org/stable/examples/analysis/README.html#hydrogen-bond-analysis should have some information on how to tell the hbond analysis what donors and acceptors are.

Perhaps others can chime in with more concrete help.

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/a78fa58a-9609-4aca-a16c-8cc5f6031824n%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

orbe...@gmail.com

Jose Martín-Roca

unread,

Oct 26, 2022, 4:53:57 AM10/26/22

to MDnalysis discussion

Hi Oliver,

Thank you very much for your answer

I think I almost solve the problem. I attach you my code (see below), but there is an error that I can not solve easily. I'm using "guess_bonds()" to obtain the topology of the system, and this method return a tuple, but the hbond.run() try to access to u._topology.bonds.values that not exists and return the error below.

Anybody knows if is there a way to guess the bonds/topology of my system with the proper structure to apply hbond.run() ?

Thanks for your time

ciao

José

-------------------------------------------------------------------------------------------------------------------------------------------------------------------

AttributeError Traceback (most recent call last)
/tmp/ipykernel_7860/3301902831.py in <module>
78 update_selections=False)
79
---> 80 hbonds.run(
81 start=None,
82 stop=None,

~/anaconda3/lib/python3.8/site-packages/MDAnalysis/analysis/base.py in run(self, start, stop, step, frames, verbose)
427 step=step, frames=frames)
428 logger.info("Starting preparation")
--> 429 self._prepare()
430 logger.info("Starting analysis loop over %d trajectory frames",
431 self.n_frames)

~/anaconda3/lib/python3.8/site-packages/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py in _prepare(self)
602 self._acceptors = self.u.select_atoms(self.acceptors_sel,
603 updating=self.update_selections)
--> 604 self._donors, self._hydrogens = self._get_dh_pairs()
605
606 def _single_frame(self):

~/anaconda3/lib/python3.8/site-packages/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py in _get_dh_pairs(self)
528 # This is because u.bonds also calculates properties of each bond (e.g bond length).
529 # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787
--> 530 if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0):
531 raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. '
532 'Please either: load a topology file with bond information; use the guess_bonds() '

AttributeError: 'tuple' object has no attribute 'values'

--------------------------------------------------------------------------------------------------------------------------------------------------------------------

from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis
import MDAnalysis as mda

import matplotlib.pyplot as plt
%matplotlib inline

##############################
######## Tables of values ####
##############################

mda.topology.tables.vdwradii= {'1':1.1, '2':1.52, 'H1':1.1, 'H2':1.1, 'O': 1.52,}
mda.topology.tables.masses= {'1': 1.00794, '2': 15.9994, 'H1': 1.00794, 'H2': 1.00794, 'O': 15.999}

##############################
######## Universe variable ###
##############################

u = mda.Universe("lammps.data",topology_format="LAMMPSDUMP", dt=0.005, atom_style="id type xs ys zs")

all_atoms = u.select_atoms("type 1 2")

Hidrogenos = u.select_atoms("type 1")
Oxigenos = u.select_atoms("type 2")

###############################

########### Charges #########
###############################
u.add_TopologyAttr("charges")
Hidrogenos.charges = 0.41
Oxigenos.charges = -0.82

###############################
########### Masses #########
###############################
u.add_TopologyAttr("masses")
Hidrogenos.masses = 1.00794
Oxigenos.masses = 15.9994

###############################

########### Names #########
###############################
nameH=[]
cont=2
for ii in range(0,len(Hidrogenos)):
if(cont==2):
nameH.append('H1')
cont=1
else:
nameH.append('H2')
cont=2

nameO=[]
for ii in range(0,len(Oxigenos)):
nameO.append('O')

u.add_TopologyAttr("names")
Hidrogenos.names = nameH
Oxigenos.names = nameO

###############################
########### Residues ########
###############################
u.add_TopologyAttr("resnames")

all_atoms.residues.resnames='HOH'

###############################
########### Guess Bonds #####
###############################
u._topology.bonds=mda.topology.guessers.guess_bonds(atoms=all_atoms,coords=u.trajectory.ts.positions,box=u.dimensions)

###############################
####### Compute H bonds #######
###############################

hbonds = HydrogenBondAnalysis(universe=u,
donors_sel=None,
hydrogens_sel="name H1 H2",
acceptors_sel="name HOH",
d_a_cutoff=3.0,
d_h_a_angle_cutoff=150,
update_selections=False)

hbonds.run(
start=None,
stop=None,
step=None,
verbose=True
)

plt.plot(hbonds.times, hbonds.count_by_time(), lw=2)
plt.title("Number of hydrogon bonds over time", weight="bold")
plt.xlabel("Time (ps)")
plt.ylabel(r"$N_{HB}$")
plt.show()

Reply all

Reply to author

Forward